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Abstract—This paper reviews the absorbing boundary 
conditions (ABCs) used in the past 40-45 years to solve 
open problems with such numerical techniques as the 
finite difference or the finite element methods. Numerous 
ABCs have been proposed over the years in the various 
domains of Physics. The review is limited to the local ABCs 
that have been used in numerical electromagnetics, even if 
some were initially introduced in acoustics. More 
specifically, the review addresses the ABCs used before 
1994, which were mainly analytical ABCs, the Perfectly 
Matched Layer (PML), used in the past 20 years, and the 
recently introduced Huygens ABC (HABC). The ABCs are 
briefly described, with an emphasis on their drawbacks 
and limitations. A special attention is paid to the 
absorption of the evanescent waves which is a critical 
question in the matter of ABCs. 


Index Terms— Absorbing Boundary Condition, ABC, 
PML, Numerical Method, Electromagnetics, FDTD, FEM. 


I. INTRODUCTION 


When solving the Maxwell Equations with such numerical 
techniques as the finite difference or the finite element 
methods, the computational domain is discretized with cells or 
elements which must be shorter than both the wavelengths of 
interest and the geometrical details. In consequence, the 
physical domain where the equations can be solved must be 
finite. This is inconsistent with most problems of 
electromagnetism that are unbounded in space. Two examples 
are the radiation of an antenna in free space and the interaction 
of an incident wave with a scattering structure. 

To overcome the contradictory requirements, the absorbing 
boundary conditions (ABCs) have been introduced. They 
replace or simulate the infinite space that surrounds the finite 
computational domain. They permit the solution to be an 
acceptable approximation of the solution that would be 
obtained within a really infinite domain. The accuracy of the 
approximation depends on the features of the ABC and on 
many other parameters, and there are some limitations on the 
simulation of free space with an ABC. Especially, the ABC 
cannot replace sources which are outside it. When physical 
sources are outside the ABC, they must be replaced with 


equivalent sources inside the ABC or with a boundary 
condition combined with the ABC. The ABCs reviewed in this 
document are only the time domain ABCs. However, most can 
also be expressed in frequency domain form and can then be 
used with such frequency domain methods as the finite 
elements or the parabolic equation. The review is also limited 
to the local ABCs. There exist global ABCs where the field on 
the boundary is estimated in function of the field in the whole 
computational domain, but their computational cost is high 
and they are rarely used in time domain methods. 

There is a need of ABCs only with finite methods. With 
such other methods as the integral equation methods, there is 
no need of ABC. Because the equations to be solved are 
equivalent to the Maxwell equations, the boundary conditions, 
and the initial conditions. The infinity of space is then 
implicit. For this reason, the need of ABCs appeared only 
when the computers became able to handle finite method 
resolutions. The first paper in the literature opening the way to 
such a possibility was probably the one of K.S. Yee in 1966 
[1], which presents the FDTD method. There is no ABC in 
this paper, but it is evident as reading it that the FDTD method 
is limited to finite domains, which implicitly suggests that 
something must be done to apply the method to the numerous 
problems of electromagnetism where the physical domain is 
not bounded by enforced boundary conditions. The first paper 
that reported the use of the FDTD method to solve an 
unbounded two-dimensional (2D) problem was published by 
Merewether in 1971 [2]. The ABC was the “radiating 
boundary”, later used in 3D by Holland [3]. 

After the Merewether work, numerous ABCs appeared in 
the literature. Schematically, there has been three periods in 
the development of ABCs. The first period, from 1971 to 
1994, where several ABCs were introduced in acoustics and in 
electromagnetics. Most are analytical ABCs, which means that 
the field on the outer boundary is computed by a certain 
formula in function of the field at some locations in the 
interior of the domain. The second period, from 1994 to the 
years 2000’s, followed the introduction of the Perfectly 
Matched Layer (PML) [4]. Because of its far better 
performance than that of previous ABCs, during this period 
almost all the works on ABCs were focused on the PML. In 
the third period, corresponding roughly to the past 10 years, 
other ABC ideas have been proposed in the literature. Some 
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are of interest and lead to ABCs that can challenge the PML, 
at least in some situations. 

This review is composed of three parts that correspond to 
the three periods of the ABC developments. Firstly the 
analytical ABC period, up to 1994, second the PML period, 
and third the recent period where new ideas emerged. The 
different ABCs are briefly described with a special emphasis 
on their drawbacks and limitations. In particular, the question 
of the absorption of the evanescent fields by the ABCs is 
discussed. This is a critical question that explains why some 
ABCs perform better than others. 


II. THE ABSORBING BOUNDARY CONDITIONS IN THE 1970 
AND 1980 YEARS 


The need of ABCs appeared at the beginning of the 70’s 
when the computer technology permitted finite difference 
methods to be used. In electromagnetics, the FDTD method 
[1] was used initially for the calculation of the coupling of the 
nuclear electromagnetic pulse (NEMP) to structures. This is a 
typical problem where an ABC is needed. NEMP coupling 
then played a capital role in the development of ABCs. In the 
USA, the radiating boundary ABC described in 2D in 1971 [2] 
and in 3D in 1977 [3] was probably the most used ABC in the 
early NEMP works. In France, use of FDTD also began with 
NEMP studies, in 1976. To this end, the matched layer (ML) 
ABC was developed [5, 6]. The ML was probably used as well 
by some workers in the USA since its use was also reported in 
[7]. 

Several analytical ABCs were introduced in the 70’s and 
80’s. Most for the acoustic waves, but they could be applied to 
electromagnetic waves as well. The most known are the 
Engquist-Majda one-way wave equations [8], and the Higdon 
operators [9, 10]. Among others some are the Bayliss and 
Turkel condition [11], the Reynolds condition [12], the Keys 
condition [13], and the Liao condition [14], which are 
reviewed in the following. 


II-1. A simple and intuitive ABC for traveling waves 


A common feature to most analytical conditions is the 
assumption that the waves to be absorbed are traveling waves 
propagating with the speed of light c. In such conditions, an 
intuitive and simple ABC can be derived. Assuming that the 
propagation is in x direction perpendicular to the boundary, 
traveling waves can be expressed in the form 


E(x,t)= f,(t-x/c) 


where E is the electric field, fg is a function of a single 
variable, and ¢ is the time. Deriving (2.1) with respect with x 
and ¢ shows that the waves E(x, t) also satisfy the equation 


ðE 10E 

—+-—=0 

ox cor 
From (2-1), the field on the boundary at a given time f is the 
field that was present at the location Ab from the boundary at 
time ft-At, where Ab is the distance of propagation during Af, 


that is Ab = c At. This is schematized in Fig. 2-1 and can be 
written as 


(2-1) 


(2-2) 


E(x, .t)=E(x, —Ab,t—At) (2-3) 


where x; is the x coordinate of the boundary. Using (2-3) as a 
boundary condition the absorption is perfect for the plane 
waves at normal incidence. But for any other incidence it 
produces a reflection which can be easily predicted. To do 
this, let us consider a plane wave E; striking the boundary at 
incidence 0, with the corresponding reflected wave E, (Fig. 2- 
2). Setting x = x, = 0 on the boundary, they can be written as 


E(x, y, t) = Ey eT (2-4) 


E.G. = Ee (2-5) 


where k,= @ cos@/c and k,=@ sin@/c are the wavenumbers in x 
and y directions. 


E(x», f) 


time t 





- 
r 7 
oO timetAt 
E(x,-Ab, t-At) : 
Boundary 
X= Xp 


Fig. 2-1. Field on the boundary at current time ¢ in function of the field at 
previous time t-At. 


Enforcing condition (2-3) on the boundary means that the sum 
of the two waves satisfies (2-3). This yields 


jot JAT jø(t-At)+ j k,âb) jo (t-At)—jk,Ab 
Ege  +Epe =E),e + Ey, é 
(2-6) 


Assuming in addition that Ab is small with respect to the 
wavelength, the exponentials can be replaced with their first 
order expansion, which yields 


Ey; + Ey, = Ey, 1— j@dt+ jk,Ab)+ E,, (1— j@At— jk,Ab) 
(2-7) 


The Boundary 


x=x,=0 


Fig. 2-2. A plane wave at incidence 8 on the plane boundary. 


From which, using At = Ab/c, the reflection coefficient r = 
Eo,/Eo; is, 





en sO (2-8) 
k ,+ø@lc 
or, since k, = @ cos@/c: 
_cosĝ-1 (2-9) 
cos +1 
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The reflection is shown in Fig. 2-3 (the traveling wave curve). 
As expected it vanishes at 8 = 0, but it is total at grazing 
incidence, and equals 17 % at 45°. It should be noticed that the 
reflection (2-8) can also be obtained by enforcing the sum of 
waves (2-4) and (2-5) into equation (2-2). The discrete 
equation (2-3) is consistent with the partial derivative equation 
(2-2). 

Reflection (2-9) has been derived for pure traveling waves. 
Consider now evanescent waves. It is also possible to derive 
the reflection coefficient. The plane wave solutions of the 
Maxwell Equations whose phase propagates in direction © of 
the (x, y) plane and whose magnitude is evanescent in the 
direction 8+7/2, can be written in the form (2-4), but with 
wavenumbers (see Appendix A) 


k= ® (cosh ycos 9+ jsinh ysin 8) (2-10) 
c 


k, = (cosh xsin- jsinh y cos 0) (2-11) 
` C 

where % is the evanescence coefficient. Traveling waves are 
the special case % = 0. The reflected wave remains in the form 
(2-5) because it propagates in direction m—® and its x 
coefficient is of opposite sign (see Appendix A). The 
derivation (2-6)-(2-8) then remains valid. Using (2-10) and (2- 
11) into (2-8) yields 


„= h gcos -1+ jsinh ysin 8 (2-12) 





cosh ycos@+1+ jsinh y sin 8 


which obviously reduces to (2-9) when x = 0. For any other 
value of x, the reflection does not vanish at normal incidence. 
For example, r(0) = 1/3 with cosh% = 2. And for strongly 
evanescent waves (coshy >> 1), the reflection is total (r = 1) 
for any 9. The reflection is plotted in Fig. 2-3 for coshy = 1.5, 
3, and 10. 
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Fig. 2-3. Reflection coefficient r for a traveling wave and for three evanescent 
waves of cos% = 1.5, 3, and 10. 


In summary, by assuming that the waves are pure traveling 
waves, boundary conditions can be easily found, in the 
continuous (2-2) or discrete (2-3) forms. They can well absorb 
the traveling waves, at least at incidences close to the normal 
incidence. But the evanescent waves are poorly absorbed, 
especially the strongly evanescent ones which are reflected in 
totality at any incidence angle. The analytical ABCs reviewed 


3 


in the following also assume that the outgoing waves are 
traveling waves. Despite that they are derived using more 
mathematics, in their first order versions they are actually 
equivalent to the above intuitive ABC (2-2) or (2-3). And their 
more sophisticated high order versions can be viewed as 
successive applications of the first order versions. From this, 
they also suffer from the drawback summarized in Fig. 2-3, 
i.e. they strongly reflect the evanescent waves. 


II-2. The Engquist-Majda ABC 
Engquist and Majda [8] derived a family of ABC’s from the 
second order wave equation, which reads, by assuming 
invariance in z direction: 
18E_ FE 5 3'E 
c ot? ax? dy’ 
Inserting a plane wave propagating in direction 6, in the form 
(2-4), yields the equation of dispersion 











(2-13) 


wlc? =k? +k? (2-14) 


which can be rewritten using k, = @ sin@/c 


k, =+ fi-sin?@ 
C 


Quantity sin’@ is small v.s. one, as long as @ is not too large. 
This suggests using an approximation of the square root. The 
simplest reads, by retaining the + sign 


(2-15) 


k,=olc (2-16) 
which can be rewritten as 
jk, = Jo (2-17) 
c 


Using the derivatives of waveform (2-4) shows that (2-17) is 
nothing but the equation of dispersion of the time domain 
equation (2-2), which is then the first order Engquist-Majda 
equation. Using it as a boundary condition produces the 
reflection (2-9). 


Boundary i= ib 





ib-1, j+1 


Fig. 2-4. Finite difference grid. 


Using (2-2) as an ABC in a finite method consists in 
discretizing it in order to permit E to be computed on the 
boundary in function of E at inner nodes and, or, at previous 
time steps. Consider the finite difference method (FDTD) with 
nodes separated with Ax and Ay in space (Fig. 2-4), sampled in 
time with step At. The discrete values of E are computed at 
nodes (i, j) and at times n. The discrete condition (2-3) is 
consistent with (2-2), but it cannot be applied as it stands since 
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in general Ab is not a multiple of the step Av. Interpolations 
would be needed to enforce (2-3). A discrete condition that 
can be used with any step Ax can be obtained [8] by a finite 
difference approximation of (2-2). More precisely, by 
discretizing the derivatives on space and time in (2-2) at 
location Ax/2 from the boundary, and at an intermediate time 
between n and n+1 (called the “box” discretization). This 
permits use of centred derivatives in space and time. With 
usual notations of the FDTD method, this reads: 








(Enj + Ex i= (Big + Eig t2 
Ax 
vl Ei Le eee ae Zo (2-18) 
c At 


from which the unknown E field at time n+1 on the boundary 
is found in function of known fields: 


E 7 Ein +w Eis =w Ei, j (2-19) 
where 
ya c At — Ax (2-20) 
cAt+Ax 


As expected, it can be easily verified by inserting (2-4) and (2- 
5) into (2-19) that its reflection coefficient is (2-9) for 
traveling waves, and (2-12) for evanescent waves. 

Higher order approximations of the square root in (2-15) 
yield higher order ABCs. The second order is obtained by 
using the first order Taylor approximation: 


k, =2(1-5sin?0) (2-21) 
c 2 
that can be rewritten using k,=@ sinO/c 
-@ +ack, anc Sk! (2-22) 


From the second derivatives of (2-4), this is the equation of 
dispersion of: 


E FE 1 , FPE 
—+¢— c — 
ot? oxot 2 dy? 
which is the second-order Engquist-Majda equation [8]. It 
should be noticed that (2-23) has been derived here assuming 
invariance in z direction. For a true 3D case, the other 
transverse derivative, i.e. 5°E/5z, must be added to 8°E/8y" in 
(2-23). 
The reflection coefficient of (2-23) can be obtained as 
previously by enforcing the sum of (2-4) and (2-5). This 
yields, for the traveling waves 


cosO-1)" 

r= — 

cosĝ+1 
which is just the square of (2-9). The second order equation 
(2-23) allows then a better absorption of the traveling waves. 
At 0 =45° the reflection is about 3% in place of 17%. For the 


evanescent waves, enforcing (2-4) and (2-5) with (2-10) and 
(2-11) into (2-23) yields 


(2-23) 


(2-24) 





__| cosh gcosĝ-1+ jsinh ysin 8 ? (2-25) 
cosh y cos +1+ jsinh ysin 0 


which is also the square of the first order reflection (2-12). 

The finite difference discretization of (2-23) can be found in 
[15] or in textbooks [16, 17]. A third order equation was also 
derived in [8] using a Pade approximation of the square root, 
achieving a reflection equal to the cube of (2-9). Later, 
Halpern and Trefelthen [18] derived other equations by using 
other approximations of the square root, such as least-squares 
or interpolation at Chebyshev or Neuman points. These 
approximations permit the angle of exact absorption (where r 
= 0) to be shifted to non zero angles (two angles where r = 0 
for the second order equation). This improves the absorption 
at wide angles. Discussions and experiments on the effect of 
the choice of the approximation can be found in [19, 20]. 
However, whatever may be their order and their complexity, 
all the approximations only considered the absorption of 
traveling waves. This means that evanescent waves are 
reflected, partially or in totality, depending on their 
evanescence coefficient. This is why in electromagnetics, 
where there are always evanescent fields, the absorption was 
not significantly improved by using Engquist-Majda equations 
of order higher than two. For a decade, up to 1994, the most 
used ABC has been the second order ABC (2-23), with its 
FDTD implementation in [15], known as the Mur ABC. The 
Engquist-Majda ABC has been also used in frequency domain 
form with the finite element method [60]. 


II-3. The Higdon ABC 


The theory of the Higdon operators was presented in two 
papers [9, 10]. Conversely to the Engquist-Majda ABCs, 
which are analytical approximations of the wave equation that 
are discretized for application to finite methods, the Higdon 
operators are directly defined in the discrete space. Higdon 
first defines two operators, K and Z L which shift backward 
the physical quantity, by one space step and one time step, 


respectively. In electromagnetics, this reads, for the E field: 
KE; ,=E; and Z“ E, =E77 (2-26) 


i-l, j 


where i is the index of the mesh in the direction of the space 
shift perpendicular to the boundary. He then considers 
operators composed of the addition of the identity operator 
with a polynomial combination of K and Z', in the form: 


B(K,Z')=1-@,K-@,Z''-a,KZ'-a,K° -a,Z” -..... 

(2-27) 

where K” and Z” mean that the shifts are applied m and n 

times (shifts by m Ax or n At). Applying B(K, Z') to the 
unknown field value on the boundary Ep 

B(K,Z™) E, =0 (2-28) 

gives Er in function of E at interior nodes of the grid (i < ib) 


at time n+1 and at nodes (i < ib) at previous times. 
Consider for instance the operator 


B(K,Z"')=1-KZ"'-wK+wZ" (2-29) 
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with w from (2-20). Its application to ES yields (2-19). The 


first order Engquist-Majda approximation can be viewed as a 
Higdon operator. This operator is then consistent with the 
equation (2-2). And its reflection coefficient is (2-9) for 
traveling waves, and (2-12) for evanescent waves. 

An important result of the Higdon papers is the possibility 
of multiplying the operators, where the multiplication is 
defined as with polynomials. Then, the reflection coefficient 
of the resulting operator is the product of the reflections of the 
two operators. This permits easy design of high order 
operators with desired features by combining low order 
operators. In fact, in his paper [10], Higdon considers the 
operator (2-29) with the following quantity w- 


_ cCAt — cos 0, Ax 


jie ON (2-30) 
c At + cos @, Ax 


Then, (2-29) corresponds to the discretization of the 


continuous equation 


ðE cos6, ðE 
poe a = 


0 (2-31) 
ox c oat 





The interest of the introduction of cos@g is that the reflection is 
no longer (2-9), it becomes 


_ cos 0 — cos 0p (2-32) 


cos 0 + cos 6, 


that vanishes at the angle 0 = Qp in place of 0 = 0. An operator 
of order p can be obtained by multiplying p operators (2-29) 
with different angles of zero reflection 0 
p 
B(K,Z")=[[U-KZ*-w,K+w,Z") 2-33) 
jal 
where the w; are given by (2-30) with the desired oj. The 
reflection is then the product of p reflections (2-32). This p- 
order operator is consistent with the continuous equation 


pP 6.. 
a= = 2)e-0 
ox 


a c ot 


(2-34) 





In practice of electromagnetics, the operator that was used 
was the second order one (p = 2), which reads, from (2-33): 


B,(K,Z')=1-(w, +w,)K +w,w,K? 
+(w, +w,)Z 1-2 +ww,)KZ1+(w,+w,)K?Z" 


+ww, Z” —(w, +w,)K Z? +K? Z” (2-35) 


The expression of E”*' on the boundary follows immediately 


by using (2-35) into (2-28). With the FDTD method, in 
general w; and w, were chosen such that 09, = 902 = 0, with 
then reflection coefficient (2-24)-(2-25). 

What should be noticed is that the second order Higdon 
ABC only needs field values at nodes located on the line 
normal to the boundary. Conversely to the discretization of 
Engquist-Majda (2-23) where values at nodes out of this line 
are needed due to the second order derivative in direction 
parallel to the boundary. This is an advantage of Higdon ABC 


in the edges of computational domains where (2-23) cannot be 
used and must be replaced with the first order (2-2). 

Since the theoretical reflection from the Higdon operator is 
the same as that from the Engquist-Majda ABC, at least when 
all the zero reflection angles 0;equal zero, the performances of 
the two ABCs are quite similar, which means good absorption 
of traveling waves, but strong reflection of evanescent fields 
in both cases. 


II-4. Other anlytical ABCs 


The first ABC reported in the literature is probably the 
radiating boundary [2, 3]. It consists in estimating the field on 
the boundary in function of the field at inner nodes assuming it 
varies as 


E(R, == (2-36) 


(t—R/c) 
R 
where R is the distance from a point somewhere in the center 
region of the computational domain. Obviously, this field 
form comes from the behavior of the field radiated far away 
from such sources as antennas or scattering structures. The 
form (2-36) is close to (2-1) with in addition the coefficient 
1/R to take account that the field decreases with distance. 
Also, the estimate relies on field values located on lines from 
the current node on the boundary to the center of the domain, 
instead of on lines perpendicular to the boundary. In principle, 
this could be more accurate than the first order Mur or Higdon 
ABCs, because of the better match of the physical field by the 
assumed evolution (2-36) in many problems. However, the 
implementation was complicated, especially in 3D [3]. And it 
seems that this ABC suffered from late time instability [21]. 
And, as Mur or Higdon ABCs, the absorption of evanescent 
fields was not addressed, which probably explains why the 
ABC had to be placed some distance from the scattering 

structure, as reported in [3]. 

Another ABC used in the early times of the FDTD method 
was that used by Taflove in his pioneering work on 
bioelectromagnetism [22, 23]. On the boundary he enforced 
the average of the fields that were present at the closest 
interior three nodes two time steps before. If only the closest 
node were used, since he worked with 2 c At = Ax the 
condition would be exact for plane waves at normal incidence, 
like with the ABC schematized in Fig. 2-1. The averaging 
removes the exactness at normal incidence, but probably 
improves the absorption at wide angles. 

The Bayliss and Turkel condition [11] is an ABC 
introduced in 1980 for electromagnetics. It can be viewed as 
an improved Merewether ABC. It also relies on the decrease 
of the field with distance from a radiating object, but it uses 
more sophisticated waveforms. Instead of (2-36) it uses the 
expansion of the field with a series of terms decreasing as the 
powers of 1/R. Operators that annihilate the outgoing field are 
derived, with different orders of approximation. The method is 
in principle more accurate than that of Merewether, but it is 
also. expressed in spherical coordinates, so that its 
implementation in a FDTD Cartesian grid is not natural and 
complicated. Its use on a circular boundary in the frequency 
domain FEM is reported in [60]. As Enguist-Majda and 
Higdon ABCs, the Bayliss and Turkel ABC only addresses 
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traveling field and then cannot work properly in the vicinity of 
sources where evanescent fields are present. 

The Reynolds ABC [12], published in 1978, is interesting 
because it connects the Engquist-Majda ABC, published in 
1977, with the future Higdon ABC, published in 1986. In 
facts, the operator (2-29), generally known as the Higdon 
operator, was previously used in the work of Reynolds. By 
reasoning as in the above section 2.1 he derived the first order 
equation (2-2) and its reflection coefficient (2-9). Then 
through a relatively complex derivation where a square root is 
approximated, he arrived at the following approximation of 
the wave equation 


19E PE 18E 
aS Fie ay aa 
coxot ox” 2 dy 

He could not obtain stable finite difference approximation to 


(2-37) and then replaced the second order derivative in the 
transverse direction y using (2-13). This gives: 


18E 28E FE 


0 (2-37) 




















bE 4 = (2-38) 
e dt K c dxdt $ ax? 
which can be rewritten as 
(2442) (2442) E=0 (2-39) 
ox cot) \ox cot 


We firstly notice that (2-39) is nothing but the second order 
Higdon equation, i.e. (2-34) with p=2 and 0)=0. Secondly, 
using the wave equation (2-13) to eliminate 5°E/8°x, equation 
(2-38) becomes (2-23). All this shows that the second order 
Engquist-Majda equation, the second order Reynolds 
equation, and the second order Higdon operator, are consistent 
and correspond to the same approximation of the wave 
equation. This can also be viewed by deriving (2-37) or (2-38) 
from (2-22) which is the equation obtained by approximating 
the square root in (2-15). For (2-38) for example, (2-22) can 
be rewritten, by replacing k, with kx: 


(2-40) 


This is the equation of dispersion of (2-38), which is then 
obtained with the same approximation of (2-15) as that used in 
the derivation of the Engquist-Majda equation (2-23). And 
consequently the Higdon operator (2-39) also can be viewed 
as derived from the same approximation. Reading [12], it can 
be noticed that Reynolds introduced a parameter p in (2-37) 
which is then present in one parenthesis of (2-39) where 1/c is 
replaced with p/c. In fact, this parameter is equivalent to the 
cos®y in the Higdon equation (2-31). It cancels the reflection 
for cos® = p, as can be verified in the table in [12]. 

The Keys absorbing equations [13] are another 
approximation of the wave equation and can be derived in a 
similar way as the Engquist-Majda equations, by starting from 
(2-14). To this end, let us rewrite (2-14) as 


@lc— ki +k; =0 (2-41) 
and let us assume that : 
Jk? +K? =Ak, + uk, (2-42) 


This gives 


@lc-Ak,- uk, =0 (2-43) 


Using the derivatives of (2-4), this is the equation of 
dispersion of 


OE 2E,13E_ 2-44 

ie Oe ot =n om 
which differs from the first order Engquist-Majda equation (2- 
2) by the presence of derivatives in both x and y directions. 
Using k,=@ cos@/c and k, = @ sin®/c shows that (2-42) is true 
only for the two angles O, and 0, that satisfy 1 cos® + u cos® = 
1. For these angles, (2-44) is exact. By expressing A and u in 


function of O; and 65, it is possible to arrive at the first order 
equation of Keys [13], which reads 











cos@, ge +sin 8, ap pe BoE =0 (2-45) 
ox oy c oat 
where 
patih y gahh (2-46) 
” 2 


By using (2-45) as a boundary equation, the reflection 
vanishes for incidences 9, and @, where (2-45) is exact. The 
two angles can be chosen at will using (2-46). A higher order 
equation can also be derived in the same way [13], starting 
from the following approximation of (2-41): 


lwie~Ak, -mk, (le Aak, = tk, }=0 


There is no detail in Keys paper on how (2-45) can be 
implemented in a finite difference method. And apparently the 
Keys equations have never been used in electromagnetics. 

The Liao extrapolation ABC [14] is a condition which relies 
on another principle. It does not assume that the waves are 
traveling waves propagating with the speed of light. Instead, it 
only relies on extrapolation of the field. For each node of the 
boundary where the field has to be predicted, a Newton 
backward-difference polynomial is found from values of the 
field at a set of nodes located on the line normal to the 
boundary, at previous time steps. Usually, such a polynomial 
is used to interpolate between values used to find its 
coefficients. Here it is used to extrapolate the field on the 
boundary at the desired future time. This results in an ABC 
which is claimed as more effective than the second order 
Higdon or Mur ABCs. This may be due to the fact that it is not 
assumed that the outgoing waves are purely traveling waves. 
From this, the absorption of evanescent waves may be better 
than with Mur or Higdon ABCs. There is no sufficient 
experiment in the literature to clearly conclude on this issue. 
The Liao extrapolation suffers from instability, but it seems an 
effective remedy was found consisting in the introduction of 
small losses. 

Strictly speaking, the Mei-Fang superabsorption [24] is not 
a new ABC, but rather a technique to improve the 
effectiveness of existing analytical ABCs with the FDTD 
method. It consists in computing the future E field on the 
boundary and the future H field half a cell from it by using the 
ABC as if the H node was on the boundary. Assuming then 
that the outgoing wave is plane and at normal incidence 
permits a relationship to be found between the errors of the 


(2-47) 
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two estimates, from which the error produced by the estimate 
of E on the advance of H with the regular FDTD equation can 
be eliminated. The H field half a cell from the boundary plays 
then the role of an exact ABC, for the assumed incident plane 
wave. In practice, it has been reported a_ significant 
improvement to the absorption of the first order and second 
order Mur ABCs. However, the improvement may not be 
better than the one that could be obtained by just increasing 
the order of the original ABC. And the technique relies again 
on the assumption that the incident waves are traveling waves, 
which means that it has probably no effect on evanescent 
fields. 

An interesting ABC that was briefly reported in the 
literature is the Betz-Mittra ABC [25]. This is a modified 
Higdon operator which is designed for absorbing the 
evanescent fields. To this end, the equation (2-2) which is 
consistent with the first order operator (2-29) is replaced with 


3E 10E 


SE pa =0 (2-48) 
ox 5 c Ot TA 


whose discretization gives the following operator 


B(K,Z™)=I-w, KZ“ -w, K +w, Z` (2-49) 


where: 


_ cAt+Ax—a@AxAt/2 . 


- _ cAt—Ax-@AxAt/2 . 
|! cAt+Ax+Q&AxAt/2 


w, = 
cAt+Ax+a@AxAt/2 


_ cAt—Ax+a@AxAt/2 
~ cAt+Ax+a@AxAt/2 


The coefficient @ in (2-48) permits the spatial derivative to 
have a contribution that takes account of the decrease of the 
magnitude of the field in the direction normal to the boundary, 
that is to take account of evanescent waves. By enforcing the 
sum of (2-4) and (2-5) with wavenumbers (2-10) and (2-11) 
into (2-48), the following reflection is obtained 


(2-50) 


3 


p2 cosh ycosð—1+ jsinh ysin + jæ/ æ 


= (2-51) 
cosh ycosð +1+ jsinh ysin- jæ/ @ 





which obviously reduces to (2-12) when &œ = 0. We can note 
that for traveling waves (cosh% = 1, sinhy = 0) r equals (2-9) 
at high frequency (@>>Q) and r = 1 at low frequency 
(@<<q). This is not a drawback since (2-49) can be 
multiplied by the normal first order operator to absorb the 
traveling waves, as in [25]. For the case where 0 = 0 and 
coshy >> 1 (strongly evanescent waves), there is no 
absorption since r = | (ratio of complex conjugates). A more 
interesting case is when 0 = - 7/2 which corresponds to the 
propagation of the phase in direction parallel to the boundary, 
and the direction of evanescence perpendicular to the 
boundary (Appendix A). In that case: 

pa ia /sinh 7 + jalo (2-52) 

1- jsinh y - jæ / @ 

For strongly evanescent waves, i.e. sinh% >>1, the reflection is 
null on condition that & =@ sinhy. It will be seen further in 
this document that a similar condition must hold as well for 
cancelling the reflection of evanescent waves by a CFS-PML. 


The condition is approximately true for the spectrum of the 
evanescent waves present in many physical problems [26, 27]. 
The relation & =œ sinh% could then be exploited to set the 
optimum Q in the same way as the & parameter of the CFS- 
PML is optimized [26, 27]. 


II-5. The complementary operators 


The complementary operator ABC was introduced in 1995 
by Ramahi [28, 29]. Strictly speaking, this is not a new 
analytical ABC, but rather a method that permits the reflection 
from existing ABCs to be dramatically reduced. Its basic 
principle relies on performing two calculations with 
complementary ABCs producing opposite reflections. If one 
calculation produces the field E in the computational domain 
and the other the field Æ“, including opposite spurious 
reflections from the outer boundary, the reflection-free (say 
exact) solution can be obtained by averaging the two results: 


E+E‘ 
2 





Ee = (2-53) 


However, there are challenging questions when applying this 
idea. Firstly, finding complementary ABCs is not trivial, and 
secondly the multiple reflections of the field from the different 
outer boundaries of the computational domain and from the 
objects inside it render the exactness of (2-53) questionable. 

The problem of the complementary ABCs has been solved 
by Ramahi in an elegant manner. He showed that 
complementary operators can be designed, based on the 
Higdon operators. The complementary operators of Ramahi 
produce opposite reflections in the continuous space as well as 
in the discrete FDTD space, rigorously. In continuous space, 
the Ramahi pair of complementary equations that approximate 
the wave equation on the boundary read 





ə Il Of COS) 5 6 (2-54) 
Ox ja ox c ot 

0 &/ 3 cos; ð 

— mee J/— |F=0 (2-55) 
ar ie ) 


a c at 
They just differ from the Higdon equation (2-34) with the 
multiplication by the derivatives 6/dx and 6/dt. These 
additional terms permit the reflections to be opposite, 
rigorously. From the Higdon theory, the reflection is the 
product of the normal operator with the reflection of the 
additional operators equivalent to the additional derivatives. 


Consider for instance the equation 6d£/dx=0. Its box 
discretization, as done with (2-2) in (2-18), yields 
a = Epa + Epi = Ep; (2-56) 
and similarly, the discretization of 5E/5r=0 gives 
Ep; Fi Epa = Epi + Es i (2-57) 


In terms of operators, the two complementary ABCs can then 
be written as 


BP(K,Z")=(I- KZ" -K+Z")B,(K,Z") (2-58) 
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Fig. 1. Multiple reflections due to terminal boundaries. comer regions, and 
the scatterer (for simplicity. the scatterer is assumed to reflect the field by 
a factor u). 
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The four different permutations of boundary operators needed to 
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solution #4 


Fig. 2-5. The various multiple reflections that occur in the computational domain (left-hand part), and the solution proposed by Ramahi to cancel corner 


reflections in 2D (right-hand part). Source: Ramahi JEEE Trans. Antennas and Propagation, vol. 46, pp. 


BP(K,Z")=(I- KZ" 4+K-Z") B,(K,Z") (2-59) 
where B,(K,Z h is operator (2-33). Using the sum of the 
incident and reflected waves (2-4) and (2-5) into (2-56) yields 
the reflection coefficient of the discrete version of the operator 
ō/ðx as: 

-jæ^t 


e zeit“ (2-60) 


-jot 
e J 


1- ei~ + e fen — elk 





n=- J-e IEA gp eN L iA 
Similarly, (2-4) and (2-5) into (2-57) give the reflection of the 
discrete version of ò/ôt 
lt+ei** —¢ gem 


r =-= : =-e 
t —jk Ax —j@dt -jk.Ax ~—ja 
Lee OM I QM 


-jot 
= jk Ax 


(2-61) 





The reflections are of opposite signs, as desired, even in the 
FDTD grid. The operators 6/dx and 6/6r are complementary. 
They could be used alone to produce opposite reflections of 
magnitude one. Obviously, one had better to multiply them 
with another operator that absorbs the waves, as in (2-58) and 
(2-59) where they are associated with a p-order Higdon 
operator. Ramahi in his papers produced results with p = 2 and 
p=3. 

It should be noticed that with wavenumbers (2-10) and (2- 
11) of evanescent waves, reflections (2-60) and (2-61) remain 
valid, i.e. the two operators yield opposite reflections. From 
this, the operators (2-56) and (2-57) yield reflections -1 and +1 
for the evanescent waves. The cancellation effect of (2-53) is 
then still valid for the evanescent fields. This is a great feature 


1475-1482, 1998 © 1998 IEEE. 


of the complementary method which permits effective 
absorption of evanescent fields, as illustrated by several 
experiments in the papers of Ramahi. 

Despite the perfectness of the pair of operators, the basic 
idea expressed by (2-53) cannot be achieved in a perfect 
manner because of the multiple reflections in the 
computational domain. This is the case in the corners where 
two negative reflections produce a positive one, resulting in 
addition of the two contributions instead of the desired 
cancellation. The non cancellation also occurs for the multiple 
reflections between the wall boundary and the object inside the 
domain, as represented in Fig. 2-5. To overcome the corner 
problem, Ramahi proposed the replacement of two calculations 
with four, with the disposal of positive and negative reflections 
(i.e. of operators on the sides of the domain) represented in the 
right-hand part of Fig. 2-5. This widely reduces the overall 
reflection, but it is costly, in 2D with the four calculations, and 
in 3D where the method can be generalized with eight 
calculations. 

To reduces the overall requirements, Ramahi then proposed 
the Concurrent complementary operator method (C-COM) 
[30]. The four, or eight, domains are limited to boundary 
regions a few cells in thickness, which are solved 
simultaneously with FDTD, while the interior domain is not 
duplicated. Each boundary region is bounded with one of the 
four (eight in 3D) dispositions of operators schematized in Fig. 
2-5. And each region is connected to the central region through 
an interface, by means of the transverse derivatives on space. 
This permits the cancellation of the reflections from the sides 


8 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


and corners of the outer domain as they enter the inner domain. 
However, these reflections are then reflected back by the 
interface towards the outer boundary where they are reflected 
another time towards the inner domain with the same sign (see 
Fig. 13 in [30]). At the end, the mechanism produces a second 
order reflection. If the reflection from the operator used on the 
boundary is R, the reflection towards the inner region is R’. 
This is a serious progress. However, in theory the same result 
could be achieved by just squaring the operator. The C-COM 
may have an advantage in the fact that it also acts on the 
evanescent waves, which is not the case with usual Higdon 
operators (but may be also the case when using the Betz-Mittra 
a coefficient). 

To conclude, the complementary operator method has been 
a major achievement in ABC technology. It permits the 
performance of existing analytical methods, mainly the Higdon 
operator ABC, to be seriously boosted. It has not been adopted 
in numerical electromagnetics mainly because it appeared 
about simultaneously with the PML ABC. Otherwise, despite 
its drawbacks, it would have been widely used, especially in 
applications where a large dynamic range is needed in the 
computational domain. And it will be seen further in this 
document that the idea of complementariness of reflections 
also is useful to interpret and understand the reflection 
observed from the PML ABC. 


II-6. The Matched Layer (ML) ABC 


The matched layer ABC [5, 6, 7] consists in using a layer of 
absorbing medium placed in front of the outer boundary of the 
computational domain, as represented in Fig. 2-6. Many 
physical media do absorb the electromagnetic waves, but all 
produce a significant reflection from the vacuum-medium 
interface, because of the mismatch of their impedance with that 
of a vacuum. From this, physical lossy media cannot be used to 
realize an effective ABC. 


-7 Matched 


j Š 
/ Region of interest `^ 


Layer 


. I 
‘. with sources a 
$ 





Fig. 2-6. The matched layer absorbing boundary condition surrounding a 
computational domain. 


An approximate solution to what can be viewed as 
contradictory requirements has been found by using an 
artificial medium with a magnetic conductivity in addition to 
the electric conductivity [5, 6, 7]. Provided a certain 
relationship holds, the impedance of the medium equals that of 
a vacuum. The reflection from the interface is then strongly 
reduced. More precisely, it equals zero at normal incidence and 
remains small as long as the incidence is not too large. 


The matched layer medium is therefore defined as a medium 
where the equations governing E and H fields read 


RL ee ie (2-62) 
ot 
p28 aetna (2-63) 


ot 


It can be easily found [6] that its impedance Z = H/E equals 
that of a vacuum, provided that 


o_o 
Eo Lo 


(2-64) 


With such a matching condition, the absorption in the medium 
over range x is 


——x 
Az=e ** (2-65) 

and the reflection coefficients r, and r, of a traveling plane 

wave at incidence 8 on the vacuum-medium interface read 


2 
isin? (1-7 2 — cos 
Ey (2-66) 
r,=—-r,=r= 5 
1-sin? IGA + cos@ 
Ey 


p Ss 
Formulae (2-65) and (2-66) are derived in [6], but they can 
also be obtained as special cases of PML media from [4] or 
[26]. At high frequency, when 6/€) @ << 1, the reflection (2- 





66) vanishes for any incidence ©. Conversely, at low 
frequency, when 6/€) @ >> 1, it becomes 
_!-cos@ (2-67) 
1+cosé 


Unfortunately, using the medium for the absorption of 
outgoing waves in numerical methods corresponds to the low 
frequency case. This is because a sufficient attenuation must be 
achieved in the ML over a thickness of a few FDTD cells or a 
few FEM elements. That is over a thickness shorter than the 
wavelength A, which is large v.s. the discretization on space in 
numerical methods. From the attenuation (2-65), this means 
that © must be large enough in order that A o/€) c > 1, or 
equivalently, o/€9 @ > 2 a. The reflection (2-66) from the 
interface of the ML used as an ABC is then, in practical 
situations, given by (2-67). This reflection is the same as that 
from a first-order analytical ABC (2-9). In addition to (2-67), 
there is the reflection from the outer PEC, which from (2-65) 
equals exp(-2 6 6 /€p c) at normal incidence for a ML thickness 
5. But it can be reduced at will by means of o and ò. For the 
evanescent waves, the reflection from PML media is discussed 
further in this document. The discussion also applies to the ML 
case (especially at normal incidence where ML and PML are 
identical). As the PML, the ML may strongly reflect the 
evanescent fields. 

The ML ABC has been used extensively in NEMP 
calculations, in France and may be also elsewhere [7]. It was 
implemented in the FDTD computer codes as a layer of only 
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three or even two cells in thickness (the computational 
domains were typically < 100° cells in size in the 70-80’s so 
that the room for the ABC was limited). The conductivity o 
was growing from a small value in the interface to a large 
enough value on the outer side, in order that the theoretical 
attenuation (2-65) over two ML thicknesses was about 99% 
(i.e. the reflection from the outer PEC was about 1%). 
Comparisons were later performed with the Higdon operators 
in this kind of applications. They showed similar 
performances, even with the second-order Higdon operator. In 
both cases the results of calculations were correct on condition 
that the ABC is placed at a sufficient distance from the 
scattering objects. This is because both the operators and the 
ML strongly reflect the evanescent waves. Some comparisons 
of the ML ABC with the Higdon and Mur ABCs can be found 
in papers [4, 31, 33]. 
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Fig. 2-7. E field produced on the PEC surface of a 64-79-20-cell aircraft 
structure with a 2-order Higdon ABC placed several distances from the 
aircraft, for a double-exponential 200-ns pulse and for a unit-step pulse. For D 
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II-7. Summary on the analytical ABCs 


There have been numerous absorbing boundary conditions 
introduced in the 70’s and 80’s. They are effective for the 
absorption of traveling waves, but little effort was paid to the 
absorption of evanescent waves, so that most ABCs reflect in 
totality the evanescent fields, except the complementary 
operators and the Betz-Mittra ABC. In many applications the 
consequence was the need of placing the ABC some distance 
away from the region of interest, in order that the evanescent 
fields decrease naturally to a sufficiently small level. This is 
illustrated with the example in Fig. 2-7 from [31] where the 
field on a airplane structure has been computed for two 
incident plane waves, with the Higdon second order ABC 
placed several distances (10, 20, 40 FDTD cells) from the 
object (from the parallelepiped that just fits the airplane). A 


distance larger than half the size of the object is needed in 
order that the low frequency evanescent fields, reflected in 
totality by the ABC, decrease naturally to a negligible level. 
The same observation can be seen also in [3] with the radiating 
boundary ABC. 


IJ. THE PERFECTLY MATCHED LAYER 


The perfectly matched layer ABC was published in 1994 [4] 
as an ABC for the solution of the Maxwell equations with the 
FDTD method. It was rapidly extended to other numerical 
techniques, and to other media than vacuum. In the following, 
section 3-1 recalls how the PML was elaborated in the simple 
two dimensional case. Section 3-2 briefly reviews the different 
forms of 3D PMLs which are well known and well 
documented in the literature and in textbooks. The other 
sections discuss with more details the critical question of the 
absorption and reflection of the evanescent waves by PMLs, 
especially in the discrete space of the finite methods. 


IIJ-1. The PML in two dimensions 


The matched layer was simple to implement, with no 
problem of stability, and was almost as effective as high-order 
analytical ABCs. At normal incidence it was perfect, so that it 
appeared that a challenging question was how to reduce its 
reflection at oblique incidence. This was achieved by 
observing that some modifications of the equations permit the 
free propagation of waves in the direction parallel to the 
vacuum-medium interface, i.e. at grazing incidence. Let us 
consider the 2D z-invariant case where the field components 
are E,, Ey, H,. The ML equations read 














n roe G-1a) 
ot “dy 
OE oH 
E) —-+ OE, =-— (3-1b) 
Cope ae 
JE 
ce eee aes (3-1c) 
ot “dy ax 


where 6 and o* satisfy (2-64). Let us now consider the 
following modified equations 


dE, (H +H) 











The (3-2a) 
A ay 
; O(H +H, 
e +0 p, = etia) (3-2b) 
0 5 dE, 
wpap g (3-2c) 
a ee 
oH 
fie eee (3-2d) 
ot oy 


where the H, component is split into two subcomponents H,, 
and H,, and where equation (3-1c) is split into (3-2c) and (3- 
2d). The medium (3-2) is represented in Fig. 3-1 with two 
plane waves propagating in the directions normal to the 
interface and parallel to the interface, respectively. Since the 
incident wave A is invariant in y direction, its propagation in 
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the medium is only governed by (3-2b) and (3-2c). For wave 
A, the medium is a matched medium, with zero reflection from 
the interface. Consider now plane wave B, assumed as set as an 
initial condition invariant in x direction, i.e. with derivatives in 
x direction equal to zero in both media. The propagation is 
governed by (3-2a) and (3-2d). The other two equations play 
no role. For wave B, the medium is viewed as a vacuum, in 
other words, at grazing incidence the medium (3-2) produces 
no reflection. Then, r = O at both normal and grazing 
incidences. 

Derivations in [4] show that r(9) = 0 at any incidence 8. The 
medium (3-2) is perfectly matched to a vacuum in the sense 
that it produces no reflection of incident waves. 


$ Vacuum Modified Matched Layer 


— oc oe 
Wave A dues Pear 











aE 
ayo H aad (3-3c) 

Bia ROR Tae 

ee erty eee (3-3d) 
ar yoy dy 


where (Ox, 6,*) and (6,, 6,*) satisfy the matching condition (2- 
64). Equations (3-2) are just the special case of (3-3) with 
0,=0,*=0. Similarly, an interface perpendicular to y direction 
produces no reflection provided that 6,=0,*=0. Combination 
of the two special cases permits an ABC to be constructed on 
the outer boundary of a 2D computational domain, as 
represented in Fig. 3-2. In the corners, both (0,, 6,*) and (0), 
0,*) are present, in such a way that the transverse 
conductivities are equal at all the small corner interfaces. This 
permits the reflection from these interfaces to be null [4]. 

An important interpretation of the PML medium (3-3) 
comes by rewriting it in frequency domain, that is by replacing 
the derivatives on time with j œ. Then, (3-3c) and (3-3d) can be 
added so that set (3-3) reduces to three equations: 











t. ieee (3-4a) 
i, Sy dy 
. | J d J a J d 0 
seeeees. e e e ee ye a seeeeeees 1 H 
; we OMY (3-4b) 
x Wave B JOEL, s, ox 
Fig. 3-1. Plane waves at normal and grazing incidences with respect to the 1 OE 1 3E 
interface between a vacuum and the modified matched medium (3-2). jo Llo H.=- + = (3-4c) 
© s ox s Oy 
More generally, the 2D (TE case) PML medium is defined Where : 
i O, 8 o 2 
dE OH. +H.) s,=14+—* > s'=14+—"_ (x,y) G4) 
Ey = +6, a (3-3a) 5 joe, jOE 
7 Equations (3-4) are identical to the Maxwell equations in a 
dE, o(H +H.) 3-3b vacuum. They just differ with the coefficients s,, 5,* and sy, s,* 
Eo Jr +O, Re aes (3-3b) in front of the derivatives in x and y directions, respectively. 
The PML can be viewed as a stretch of coordinates by the 
PML (6415010952) PML(0,0,0,5,0,5) PML (0,200,203) 
Perfect Conductor a Vacuum „= 
PML(,,.0,,,0,0) PML(©,2,0,2,0,0) 
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Fig. 3-2. The 2D PML ABC composed of PML media ensuring no reflection at all the interfaces. 
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complex factors s,, and s,,* which are equal when the matching 
condition (2-64) holds. For zero reflection at an interface, the 
medium is only stretched in the direction perpendicular to the 
interface, that is x with (3-2), and left unchanged in the other 
direction. This is the fundamental difference with the matched 
medium (2-62) which can also be viewed as a stretched 
medium, but in all the directions of space. The interpretation of 
the PML in terms of stretched coordinates is due to Chew and 
Weddon [32]. It is essential in view of extending the PML to 
other media than vacuum. 

Using (3-3), the equation of dispersion in the PML can be 
easily derived by inserting a waveform like (2-4) 


ok 


k? 
Fe 





= - (3-5) 
C SSe SS 
It is satisfied by the following wavenumbers 
k, = Js 5° cosd (3-6a) 
c 
k, = [s,s sinð (3-6b) 
P a NEE 
ô 
rr 
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Fig. 3-3. The reflection of a plane wave from the wall PML. 


In medium (3-3) only stretched in x direction (s,= s,*= 1) the 
wavenumber k, is the same as in a vacuum. From this, the 
equation k,,=k,. that holds at the vacuum-PML interface 
implies that the angle of propagation of an incident wave is left 
unchanged when it penetrates into the PML (in other words, 
the Snell law reads 0, = 05 at the vacuum-PML interface). 
Using then k, (3-6a) into (3-3) shows that the medium remains 
an absorbing medium. More precisely, in the PML where 6,= 
0,*=0 and (6,, 0,*) satisfy (6-24), the field E (or H) at distance 
x from the vacuum-PML interface can be written as 


O, cos 0 
ee 
Eo c 


E e (3-7) 


vacuum 


PML 7 


which means that the waveform in the PML is the same as in a 
vacuum, with just in addition the attenuation by an exponential 
factor. We notice the attenuation depends on the angle 0. At 
normal incidence it equals the attenuation of the matched 
medium (2-62). It vanishes at grazing incidence, which is not a 
drawback for using the medium as an ABC. From (3-7), the 
reflection produced by the side PMLs (also known as the wall 


PMLs) of thickness 6 (Fig. 3-3) reads, in the case where the 
conductivity varies with x: 


cos 8 





Í (x)d: 
0, (x) dx 
Ec 10 * 


Rce (3-8) 


In principle, R(®) can be chosen as small as needed, by 
increasing ©, or ð, or both. 


HI-2. The PML in three dimensions and for any inner media 


The above 2D PML can be extended in a trivial manner to 
3D, by splitting the six components of the field into twelve 
subcomponents, and by splitting the six components of the 
Maxwell equations into twelve sub-equations [33]. This 3D 
counterpart of the 2D case (3-3) is known as the split PML 
where there are then twelve quantities (the sub-components) to 
be advanced in time (only ten in the wall PMLs). Another 
method [32] consists in generalizing (3-4) to 3D, that is 
stretching the coordinates in the Maxwell equations by 
coefficients like (3-4d). From this, other time domain 
equations can be derived, namely the convolutional PML 
(CPML), and the near PML (NPML). In the CPML introduced 
in [36] there are only the six components of the field, but 
twelve additional quantities (four in wall PMLs) must be 
computed and advanced in time by means of auxiliary 
equations. The NPML is an approximation to the PML which 
is identical to a PML only when the conductivity is uniform in 
the medium |37], but it also produces no reflection even when 
the conductivity varies [38, 39]. It also involves twelve 
additional quantities governed by twelve auxiliary equations. 

Another absorbing medium which produces no reflection at 
the interface with the inner medium has been proposed [34, 
35]. It is known as the uniaxial PML. There the fields E and H 
are stretched with the coefficients (3-4d), in place of the 
coordinates. This also results in zero reflection [34]. Its time 
domain equations are those of an anisotropic medium with six 
additional quantities (two in wall PMLs) and six auxiliary 
equations [35]. This medium can be used easily with the finite- 
element method. 

The PML medium can be used as an ABC on the outer 
boundary of a computational domain, by extending the 2D case 
in Fig. 3-2. By an adequate choice of the conductivities [26, 
33] there is no reflection at all the interfaces in the problem 
space, i.e. at the vacuum-PML interfaces as well as at the inner 
interfaces between the different PMLs in the corners and 
edges, where transverse conductivities always are identical on 
the two sides of the interfaces. 

The PML can be extended to any media by applying the 
stretched mesh method (p. 38-42 in [26]), especially to lossy 
media, anisotropic media, or non homogeneous media. And 
more recently, to metamaterial media. 

Each time domain PML has its own advantages and 
drawbacks, in terms of needed storage, computational time, 
and easiness to implement with complicated media. A 
comparison of the computational requirements can be found in 
(p. 68, [26]). From the works reported in the literature, it seems 
that the most popular implementation in the FDTD method is 
the CPML which can be more easily extended to any media 
than the split PML (this is also true for the NPML). For 
detailed descriptions of the different PML media and of their 


12 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


time domain implementations in the FDTD method, the reader 
is referred to the existing textbooks [16], [17], and [26]. The 
PML ABC can also be used with such frequency domain 
methods as the parabolic equation method [61] and the finite 
element method where its most used version is the uniaxial 
PML [34]. 


IIJ-3. The evanescent waves in the PML 


In the above and in most papers in the literature, especially 
in [4, 32], only traveling waves absorbed with the coefficient 
in (3-7) have been addressed. They correspond to the 
wavenumbers (3-6). However, it can be easily verified that the 
following wavenumbers are also solutions of the 2D equation 
of dispersion (3-5): 


k= eS s. (cosh ycos6 + jsinh ysin 8) (3-9a) 
c 


SO lee inĝ- jsi 3-9b 
k, VSS (cosh ysin @ — jsinh ycos@) (3-9b) 
where x and @ are free parameters. They just differ from the 
wavenumbers in a vacuum (2-10)-(2-11) with the square root 
of the stretching factors. And when y=0 they reduce to those of 
a traveling wave propagating in direction 9 (3-6). Inserting (3- 
9) into the waveform (2-4), and using the change of 
coordinates defined in appendix A, the following waveform is 
obtained: 


Oy 








cosh O, 5 $ 
a Zy- *— sinh ysin ox] -2 sinh yY = cosh y cos 8 x 
(a e 


jfi s 
= c Eoc 0 
E=E,e 


(3-10) 


The two terms depending on X and Y are nothing but the 
waveform of an evanescent wave in a vacuum (see appendix 
A), whose phase propagates in direction X forming the angle 0 
with x, and whose magnitude decreases in direction Y 
perpendicular to X. The two additional exponentials are a 
phase term and an absorbing term. This shows that in a PML 
the evanescent waves can be written as 


©, cosh y cos 0 
OO A ORS 
Eo c 


E (3-11) 


vacuum 








|E puz | = 


Since coshy > 1, the absorption of evanescent waves in the 
PML is larger than that of traveling waves (3-7). 

It was derived in [40] that in continuous spaces the 
evanescent waves are not reflected from a vacuum-PML 
interface. As the traveling waves they can penetrate in totality 
into the PML. The angle of the propagation of the phase © is 
left unchanged through the interface, as well as the 
evanescence coefficient %. In theory, the PML is then perfectly 
matched to the evanescent waves. However, in discrete spaces 
we can expect some spurious effects when the PML 
conductivity has been chosen in order that the absorption (3-7) 
is sufficient to absorb the traveling waves. Then, the absorption 
of strongly evanescent waves, when coshy >> 1, may be quite 
high. In consequence, it may correspond to a total absorption 
over a distance shorter than one cell with the FDTD method or 
one element with the FEM, as represented in Fig. 3-4. From (3- 
11), this occurs when: 


e 


cosh y >> 





Eoc (3-12) 
„Ax 

It is obvious that the FDTD mesh cannot properly sample 
such waves that decrease too rapidly. From this, a significant 
reflection of strongly evanescent wave can be expected from 
the vacuum-PML interface. This has been confirmed by 
numerical experiments and by the theory of the evanescent 
waves in the discrete FDTD method [40], as discussed with 
details in the next section. We can also notice that when (3-12) 
holds, the phase term in (3-10) produces rapid variations of the 
phase with respect with the cell Ax. This cannot be sampled by 
the mesh as well, and also may contribute to the spurious 
reflection. 


H-4. The numerical reflection from the PML 


In theory, any desired small reflection (3-8) can be achieved 
by using a large enough product ©, 5. More specifically, any 
small reflection can be achieved with an arbitrary short 
thickness 6, provided that ©, is large enough. This means that 
the PML could be one cell in thickness with the FDTD method 
(one element with the FEM). In practice, as soon as the early 
experiments [4], it appeared that the PML can achieve far 
better performance that previous ABCs, but on condition that it 
is some FDTD cells in thickness. There exists a spurious 
reflection, not predicted by (3-8), which is due to the 
discretization of the medium. This reflection decreases when 
the thickness grows and when 6, grows from a small value in 
the interface to a larger value in the outer cell of the PML. This 
is in accordance with the known fact that at an interface 
between two physical media, as air and a dielectric, a spurious 
reflection is added to the physical reflection when using 
discrete methods. The same occurs with the air-PML interface. 


Wave magnitude 


Homogeneous traveling wave 


X2>xXi> 0 


Strongly 
evanescent 
wave 





Fig. 3-4. Theoretical decrease of evanescent waves in the FDTD grid for 
several values of the parameter ¥. 


In the following, the discussion of the numerical reflection is 
limited to the FDTD method. The reflection with other discrete 
methods (FEM, TLM, others), would not be very different 
because it is mainly produced by the too rapid decrease of the 
field with respect with the space sampling in the PML. More 
specifically, using the same sampling (the same cell or 
element), the numerical reflections produced by different 
methods are similar, as it is illustrated with the FDTD-TLM 
comparison reproduced in [26]. 

We must now distinguish two cases, the pure traveling 
waves, and the evanescent waves, which are fundamentally 


13 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


different from the point of view of the numerical reflection 
produced by a PML. 


The numerical reflection of traveling waves 


Numerous experiments were reported in [4] and elsewhere, 
but the first theoretical study of the spurious reflection of 
traveling waves was reported in [41] for the air-PML interface. 
More generally, for a PML where the conductivity grows from 
one cell to the next, the reflections from all the inner interfaces 
must be taken into account to obtain the overall reflection. To 
predict this overall reflection, the formulae (a tri-diagonal 
system) derived in [40] can be used in the special case of 
traveling waves (coshy=1). 

In general, the numerical reflection of traveling waves is 
relatively small. It is not a challenging question. It depends on 
the profile of conductivity in the medium, and satisfactory 
results can be obtained with polynomial profiles (power 2 or 3 
usually). Fig. 3-5 shows an example. With a parabolic profile 
and a PML only four cells in thickness, reflections as low as — 
80 dB can be achieved at normal incidence. Notice that when 
the theoretical reflection R(0) decreases, the numerical 
reflection (difference between theory and FDTD in the figure) 
grows. It is rather counterproductive to use a too small R(0). 
For given PML thickness and profile of conductivity, there is 
an optimum R(0) where the actual FDTD reflection is 
minimum. Reducing further the reflection requires increasing 
the PML thickness. The experiments in [4] demonstrate the 
evolution of the actual reflection with some parameters 
(thickness, profile, R(O)). It is easy to predict this reflection 
using the formula in [40], which is also available in the form of 
software [42]. 
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Fig. 3-5. Comparison of the theoretical reflection R(®) with the reflection 


observed in FDTD experiments with plane waves. 


The numerical reflection of evanescent waves 


The reflection of evanescent waves has been, and remains in 
part, a challenging question. In the initial experiments [4], it 
appeared that in such problems as wave-structure interactions, 
strong low frequency spurious fields are present, as with 
analytical ABCs. However, by using a thick enough PML, the 
reflection could be reduced and the PML placed in the close 
vicinity of the scattering objects. Such an achievement was not 
possible with Mur of Higdon ABCs, with which the ABCs had 
to be placed far from the objects, out of the evanescent region. 
This kind of problems was initially analyzed empirically, and 
guidelines for choosing the PML were reported in [43, 44]. 


Understanding the reason of the low frequency spurious 
reflection from PMLs began with the pioneering work of De 
Merloose and Stuckly [45] who showed that the actual 
reflection from FDTD PMLs can be predicted theoretically for 
waves evanescent in the direction perpendicular to the PML 
interface (at the end of a waveguide). This reflection can be 
large and event total, because of the rapid variation of the 
phase of evanescent waves in the PML, upon characteristic 
distances far shorter than the FDTD cell size. This case 
corresponds to 8 = 7/2 in wavenumbers (3-9). 

Following the way open by [45], the theory of general 
evanescent waves in PML media was derived, both in 
continuous space and in the discretized space of the FDTD 
method [46, 40]. The theory permits the spurious reflection of 
low frequency observed in experiments to be interpreted and 
well understood. As discussed in the previous section, the 
evanescent waves penetrate without reflection in the 
continuous PML. This is no longer true in the discrete FDTD 
space. For a PML with uniform conductivity the reflection 
from the interface is given by an analytical formula (23 in [40] 
or 5-45 in [26]). In the case of strongly evanescent waves, i.e. 
when coshy, satisfies (3-12), the formula reduces to: 


: fe 
j 
paf (3-13) 
T T 
f 
with 
zo ris (3-14) 
fe 27 Ey 


where Ox is the conductivity in the interface which may be 
different from the uniform conductivity 6, set at the other 
FDTD nodes of the PML (as 0,9 = 6,/2). 
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Fig. 3-6. FDTD reflection from the interface between a vacuum and a uniform 
infinite PML, from traveling wave (coshy = 1) to strongly evanescent waves 
(coshx = 1000). 


Formula (3-13) shows that when (3-12) holds the reflection 
from the interface decreases at high frequency (f>> f.), but it is 
total (R,=-1) at low frequency (f << f,). The total reflection 
below f,, was expected since when (2-12) holds the evanescent 
waves are absorbed, in theory, over less than one FDTD cell, 
which cannot be achieved by the FDTD method. Fig. 3-6, from 
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[40] shows the reflection coefficient computed with the 
formula in [40, 26], at the incident angle 8 = 7/4, which means 
that the evanescence is also 7/4 from the interface. The PML 
conductivity in the interface is such that f. = 10 MHz. The 
reflection is plotted for cosh% varying from 1 (traveling wave) 
to 1000. For the traveling wave the reflection grows as the 
frequency is approaching the cut-off of the FDTD mesh. It is 
very small (- 100 dB) at low frequency. Conversely, for the 
strongly evanescent waves the reflection grows as the 
frequency decreases and it is total below f,, in accordance with 
(3-13). Notice that the condition (3-12) for the validity of (3- 
13) is cosh% >> 40 in the conditions of Fig. 3-6. This is well 
verified. 
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Fig. 3-7. Comparison of the spurious FDTD field observed on a 2D thin plate 
with the factor 1-R, where R is the theoretical reflection of evanescent waves 
by a PML with conductivity growing geometrically (see more details in [40]). 
The oscillations of the spurious field illustrate the reflection by the inner 
interfaces of the PML, which are alternatively of electric and magnetic kinds. 
The two parts should be compared only below 30 MHz (above, the traveling 
waves are predominant, they are absent in the theoretical lower part). 


In the case of a PML with ©, growing from a small value in 
the interface, the situation is more complex. The reflection is 
given by a simple tri-diagonal system [40, 26, 42]. For f < fo 
and strongly evanescent waves (3-12), the reflection also 
reduces to (3-13). But every inner interface has its own 
frequency f, proportional to its conductivity. So that 
frequencies that penetrate through the vacuum-PML interface 
may be reflected from the inner interfaces. Since in the FDTD 
grid there are electric and magnetic interfaces (interfaces 
where the tangential fields are either E or H, separated with 
Ax/2), the internal reflections are alternatively of opposite 
signs. From this, the numerical reflection of strongly 
evanescent waves oscillates as the frequency grows. This is 
well predicted by the theoretical formula in [40, 26], and very 
well verified in experiments. An experiment reported in [40] is 
reproduced in Fig. 3-7. The upper part shows the F field on the 
surface of a thin PEC plate stricken by an incident wave, with 


several PMLs placed close to the plate, two FDTD cells from 
it. In each case, the low frequency plateau is strongly 
erroneous, corresponding to the reflection by the vacuum-PML 
interface for f < fẹ. At higher frequencies, the oscillations are 
produced by the reflections from the successive inner 
interfaces. The lower part shows that the general shape of the 
reflection of the evanescent waves can be very well 
reconstructed by means of the formula in [40, 26]. In facts, the 
PML works as a set of complementary reflectors, in the sense 
of the complementary operators of Ramahi (see section 2-4). If 
the ratio of frequencies f, of the inner interfaces is not too 
large, each frequency experiences several partial reflections of 
opposite signs from several interfaces. The overall reflection of 
this frequency is then smaller than total reflection, and the 
closer the f,’s, the smaller the overall reflection. 

In summary, the evanescent waves are reflected by the PML 
in such discretized spaces as the FDTD space. However, and 
conversely to most analytical ABCs, it is possible to reduce the 
reflection, by increasing the PML thickness and by using a 
conductivity that grows in the PML. Two conditions must hold 
for the reflection of evanescent waves is small enough: 

1/ the conductivity Oy in the vacuum-PML interface must 
be small enough in order that the frequency f, (3-14) is smaller 
than the lowest frequency of interest of the problem. 

2/ the conductivity must grows at a small rate from one 
FDTD node to the next, in order to reduce the reflection from 
the inner interfaces, or, in other words, to maximize the 
Ramahi effect. 

Optimum PMLs satisfying 1/ and 2/ for solving wave- 
structure interaction problems were elaborated [43, 44]. 
However, they cannot be applied to any problems because they 
depend on the general form of the evanescent fields and on the 
desired accuracy. In practice, it seems that is many, not to say 
most, cases reported in the literature, the PML is chosen 
mainly empirically, by performing several tests to find the 
minimum thickness needed for solving the problem under 
interest. It should be noticed that the choice of the PML (split 
PML, uniaxial PML, CPML, NPML) has little effect on the 
FDTD reflection of both traveling and evanescent waves. This 
was demonstrated both theoretically and by experiments in 
[47] and [39]. 

The next section is devoted to the complex frequency shifted 
PML (CFS-PML), which permits the reflection of evanescent 
waves from PMLs to be dramatically reduced in many physical 
problems. This is obtained by means of a modified stretching 
factor. 


Il-5. The complex frequency shifted PML 


The CFS-PML was introduced in 1996 by Kuzuoglu and 
Mittra [48] with the intention of rendering the PML medium 
causal. To this end, they replaced the stretching factors (3-4d) 
with the following, for the stretch in x direction 


Oo 
aa 


> (3-15) 
A, + JME, 


where ©, is a real positive parameter. The stretching factor 
remains finite as frequency vanishes which render the medium 
causal [48]. Another interest of this modified stretched factor 
was found later [36, 47]. It permits the reflection of evanescent 
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waves from PMLs to be dramatically reduced in many 
applications of finite methods. 
With (3-15), it is evident that frequency 


Sa 


= (3-16) 
RR Ei 


is a critical parameter. For f >> fy, the CFS-PML is a normal 
PML since then o, is negligible in (3-15). Conversely, when f 
<< fo, the complex term is negligible and the CFS-PML is just a 
real stretch of coordinates. It does not absorb the waves. 

Using wavenumbers (3-9) and coefficient (3-15) into (2-4), 
the waveform can be obtained in the CFS-PML [47]. As 
expected, for f >> fa it reduces to the waveform in a normal 
PML, so that (3-11) remains valid. Conversely, if f << fa the 
waveform is different [47, 26]. Its magnitude decreases in the 
PML according to 


f o,sinh ysin@ | 
Ía 


E ae (3-17) 


vacuum 


e 








|E pu | = 


where 8 is negative when the evanescence is towards the PML 
(appendix A). For low frequency (f << fa) and strongly 
evanescent waves where sinhy, ~ coshy >> 1, it is evident that 
the attenuation coefficient of (3-17) is by far smaller than that 
of (3-11). This may result in a reasonable attenuation of 
evanescent waves, and then a dramatic reduction of their 
reflection from the discrete vacuum-PML interface. 
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Fig. 3-8. FDTD reflection from the interface between a vacuum and a uniform 
infinite CFS-PML, from traveling wave (coshy = 1) to strongly evanescent 
waves (coshy = 1000). 


As with the standard PML, the FDTD reflection from a 
CFS-PML with uniform conductivity is given by a formula, 
and from a CFS-PML with growing conductivity by a tri- 
diagonal system [47, 26]. In both cases, for strongly 
evanescent waves (3-12) the reflection from the vacuum-PML 
interface tends to the following formula, in place of (3-13): 

jfa 

f O vo 
j- jfa G9 + Oxo 


a (3-18) 


o0 


with 


= ao A: O o 
2A Ea 


fea 6-19) 


where 0,9 and Oxo are O, and 6, in the FDTD interface. The 
reflection (3-18) tends to G,o/(O,0+O x0) at low frequency in 
place of the total reflection of the standard PML (3-13). This is 
illustrated by Fig. 3-8 which is the counterpart of Fig. 3-6. 
With Qo = 0.01 S and the conductivity in the figure, the 
reflection of evanescent waves is about -25 dB, whatever may 
be the evanescence parameter coshy. Other results illustrating 
the FDTD reflection by the CFS-PML can be found in [47]. 

The effectiveness of the CFS-PML for reducing the 
reflection of evanescent waves is demonstrated in the 
following, for a waveguide problem and for the interaction of 
an incident wave with an object. Let us first consider a parallel 
plate waveguide with separation a between the plates. The 
mode n is evanescent below its cutoff which reads œ= n T c / 
a, and its sinh% parameter is given by 


2. 
Onoff =F 1 


sinh y =+ s (3-20) 
w 
Thus, far below the cutoff the following holds 
o sinh y = +”7E (3-21) 
a 


From this the attenuation in (3-17) does not depend on 
frequency far below the cutoff. By an adequate choice of the 
free parameter fa, the attenuation of the evanescent waves (3- 
17) and that of the traveling waves (3-7) can be set equal. 
Assuming that sin6 ~1, this is realized with the following fy: 


f= uae (3-22) 
2a 
which yields, using (3-16): 
a TCE, (3-23) 
a 
Seutoff 
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Fig. 3-9. Coincidence of the transition frequency of the CFS-PML with the 
cutoff frequency of the parallel-plate waveguide. 


We can notice that fy = femo The situation is then as 
represented in Fig. 3-9. For the traveling waves, above the 
cutoff, the CFS-PML is a standard absorbing PML. 
Conversely, for the evanescent waves, below the cutoff, the 
CFS-PML is just a real stretch of coordinates, with no 
absorption. There the evanescent waves decrease in a natural 
manner over the stretched space, resulting in an apparent 
absorption identical to that of the traveling waves, and in 
consequence with little numerical reflection from the vacuum- 
PML interface. Numerical experiments with this waveguide 
problem are reported in [49]. They show that the CFS-PML 
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can achieve reflections of evanescent frequencies below -80 
dB with a CFS-PML only 6 FDTD cells in thickness. 

Another case where the CFS-PML permits the reflection of 
low frequency to be reduced is when an incident wave strikes a 
scattering structure. In that case, a relationship like (3-21) also 
holds, more precisely (see appendix B): 


æ sinh y= pÅ (3-24) 
w 
where w is the largest size of the structure and p is a parameter 
of the order of unity. By assuming that p = 1, the attenuations 
of traveling (3-7) and evanescent (3-17) waves are equal 
provided that 
je eee (3-25) 


% 


w 


This corresponds to fo =fo/% where fo = c/2w is the fundamental 
resonance of the structure. The situation is then the same as in 
Fig. 3-9. For the high frequencies above fo, which are traveling 
frequencies, the CFS-PML is a standard absorbing PML. And 
for the low frequencies below fo, which are evanescent, the 
CFS-PML is just a real stretch of coordinates that permits them 
to decreases in a natural manner. The FDTD reflection from 
the vacuum-PML interface, and from the inner interfaces of the 
PML, is then quite small for all the frequencies. 

The CFS-PML for wave-structure interactions has been well 
studied and optimized [50, 26, 27]. It permits low reflections to 
be achieved with CFS-PMLs only 3-5 cells in thickness, with a 
short separation of two FDTD cells between the PML and the 
structure of interest. Numerical experiments and the 
recommended parameters of the CFS-PML for this kind of 
applications are provided with details in [27]. Fig. 3-10 shows 
an example, with the same aircraft structure as in Fig. 2-7, but 
re-grided with smaller FDTD cells (6.25 cm in place of 50 cm). 
The results are superimposed on the reference solution with a 
CFS-PML only 3-4 cells in thickness placed two cells from the 
structure. Achieving such a performance with the second order 
Higdon operator would require a huge amount of vacuum 
around the aircraft, as demonstrated by Fig. 2-7 where the 
results are not perfect even with a separation of the order of 
half the size of the object (this is the physical separation 
between the object and the Higdon operator which permits the 
evanescent wave to decrease, computing the equivalent of 
result D = 40 cells in Fig. 2-7 would require D = 320 cells with 
the re-grided aircraft). 

From the above, the CFS-PML permits an effective 
absorption of evanescent waves with little reflection, by 
choosing the parameter œ, in such a way that the attenuation of 
the traveling and evanescent waves are equal. However, a 
weakness of the method is that the optimum Q, may not be 
unique. With the waveguide problem, the optimum ©, depends 
on the mode (3-23). Similarly, with the scattering problem, ©, 
depends on the parameter p. The problem has been overcome 
in the optimum CFS-PML [27, 26] by using a parameter O, 
that varies in the PML, from a value larger than (3-25) in the 
interface to a value smaller than (3-25) at the outer end. This 
results in a satisfactory optimum CFS-PML for the kind of 
problems addressed in [27], mainly electromagnetic 
compatibility (EMC) problems. A similar strategy could also 
be used with waveguide problems where several modes are 
present, at least when the cutoff frequencies of the modes are 


not too different. This has never been tested. Another weakness 
of the CFS-PML optimisation is that the Physics of the 
problem must be known a priori, that is the relationship 
between frequency and the sinhy parameter of the evanescent 
waves must be roughly known. This may be difficult or not 
possible in some problems. In the literature, the CFS-PML is 
widely used. However, in most cases the choice of 0, is done 
without reasoning and without convincing argumentation. 
Probably empirically by performing some experiments with the 
problem under investigation to find which value gives the best 
results. 
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Fig. 3-10. Electric field on the surface of a 504-632-160-cell airplane with 3-, 
4-, 5-cell CFS-PMLs, for a unit-step and a 200-ns double exponential waves. 
From the interface to the outer end, the conductivity grows by a factor of 10 
and the parameter 0, decreases from 50 to 1/5 of value (3-25). The CFS-PML 
is two cells from the structures (from the end of the fuselage, the wings, the 
tail). More details can be found in [27]. 


Two remarks to conclude. Firstly, it should be noticed that 
Œ, is a physical parameter, homogeneous to a conductivity 
(same units). This strongly suggests that its value should not be 
chosen as a pure arbitrary numerical parameter, but rather in 
function of physical parameters, as done in the two problems 
addressed in the above. Secondly, the CFS-PML is in general 
implemented in FDTD as a convolutional PML (CPML) [36], 
although an auxiliary equation implementation is possible as 
well. There is sometimes some confusion in the literature about 
the origin of its interesting possibilities. The CFS-PML can 
reduce the reflection of evanescent waves because of its 
stretching factor (3-15), not because of its convolutional 
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implementation, i.e. another implementation would produce 
the same beneficial effect on the evanescent waves. 


IlI-6. Summary of the Perfectly Matched Layer 


The PML has been viewed as a great advance in the field of 
the simulation of free space. In comparison with the numerous 
existing analytical ABCs, it permits a large reduction of the 
reflection of the outgoing waves in open problems solved by 
means of finite methods. Absorption of traveling waves is 
probably not better than that of such analytical ABCs as the 
Mur ABC or the Higdon ABC which was satisfactory (see next 
section on the actual performance of analytical ABCs). The 
principal advantage of the PML is its possibility of absorbing 
the evanescent fields which are present in almost all the 
problems of electromagnetism. In facts, for the reasons 
recalled in the above, the PML ABC may strongly reflect the 
evanescent fields in the discrete space, because its theoretical 
absorption of such waves is too high. This can be overcome by 
two means. The first one is to use a relatively thick PML where 
complementary reflections (Ramahi effect) from the successive 
rows of the discrete PML result in an overall reflection of 
evanescent fields which can be as small as desired, provided 
that the PML is thick enough. The other mean is the use of the 
CFS-PML which relies on a modified stretching factor. In 
many physical problems it permits the theoretical absorption of 
the evanescent fields to be the same as that of traveling waves, 
reducing then their reflection from the discrete PML. However, 
these remedies have some drawbacks. Choosing the best and 
optimum parameters of the PML is not evident and remains 
dependant on empirical knowledge, and on a certain a priori 
knowledge of the Physics of the problem to be solved. From 
this, the question of the absorption of evanescent fields 
remains a challenging question which prevents from 
considering the PML as an ideal ABC. For these reasons, it is 
clear that further researches in the field of absorbing boundary 
conditions remain of potential interest, mainly concerning the 
question of the evanescent waves. 


IV. ABSORBING BOUNDARY CONDITIONS IN THE LAST DECADE 


Despite the fact that the PML was considered as an effective 
ABC for solving any problems of electromagnetism, some 
innovative works have been published in the literature over the 
past 10-15 years, reporting novel ideas in the field of ABCs. 
For example, the lacunae-based ABC, introduced in the 
context of acoustic waves, has been the subject of several 
papers in years 2003-2012, from [51] to [52]. This ABC 
mainly relies on the fact that the acoustic or electromagnetic 
response to a compact in time source is lasting over a finite 
duration. By partitioning any time domain source into a set of 
short elementary sources, this permits the solution of open 
problems to be computed for any arbitrary long time using a 
finite domain, without reflection from the outer boundary. The 
principle of this ABC is attractive. However, there is no 
comparison with previous ABCs in the literature papers, nor 
example of application to problems of interest in 
electromagnetism. And it seems that the overall computational 
requirements of this ABC are far larger than those of the PML 
ABC. 


In the following, we just review three ABCs published in the 
years 2000’s for the absorption of electromagnetic waves. 
More precisely, two ABCs that were published almost 
simultaneously, the multiple absorbing surface [53] and the re- 
radiating boundary condition [54, 55]. They rely on the same 
fundamental principle. And a third ABC, the Huygens ABC 
[56], which is the generalization of the previous two. The 
HABC is promising in several aspects and could challenge the 
PML ABC, at least in some applications. 


IV-1. The re-radiating boundary condition (rRBC) 


The rRBC [54, 55] was introduced for use with the FDTD 
method. Its principle consists of generating a wave opposite to 
the outgoing wave leaving the computational domain, by 
means of a Huygens surface. If the generated wave were 
exactly opposite to the outgoing wave, the cancellation would 
be perfect, as represented in Fig. 4-1. Whatever may be the 
outer boundary placed behind the Huygens surface, there 
would be no backward reflection in the computational domain. 
In the papers [54, 55], the generation of the outgoing wave 
relies on the concept of teleportation of fields. The field one 
time step in the past and one space step backwards is 
“teleported” upon the Huygens surface where it is used to 
radiate the opposite field. This outgoing field is not exact, it is 
only an estimate, so that the cancellation is not perfect. From 
this the outgoing field is partially reflected. To improve the 
absorption, the computational domain in [54, 55] is ended with 
an impedance condition equivalent to a first order operator 
ABC. 

There is neither derivation nor theoretical prediction in [54] 
and [55], only FDTD experiments are provided. The authors 
observed that the field outside the rRBC resembles the 
derivative on time of the outgoing field. They used the terms 
derivation and re-integration to characterize the direct and 
reverse effects of the rRBC. As it was derived later [56], the 
effect of the rRBC is indeed a derivation of the outgoing 
waves, and an integration of the ingoing waves. 

Several experiments reported in [54, 55] show that the rRBC 
performance can be similar to that of the PML ABC. However, 
generalization of the conclusions of these experiments is 
questionable, because they were performed in a 2D domain 
with the electric field perpendicular to the domain. It is known 
that any ABC performs very well with this kind of problems, 
because no evanescent waves have to be absorbed. Especially, 
Higdon and Engquist-Majda ABCs also yield quite good 
absorption in that case [20]. In realistic FDTD applications 
where evanescent waves are present, results with the rRBC 
would be far poorer, because the teleportation used in this 
paper produces a total reflection of evanescent waves, as 
discussed in the following. 


IV-2. The multiple absorbing surfaces (MAS) 


The absorbing surface concept in [53] is introduced in the 
FDTD method by splitting the computational domain into a 
scattered-field region and a total-field region. There is no 
explicit reference to equivalent currents radiating a field 
opposite to the outgoing field. Nevertheless, what is done in 
this paper is equivalent to the generation of a wave opposite to 
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Fig. 4-1. The principle of the Huygens absorbing boundary condition. 


the outgoing field by means of a Huygens surface, as in Fig. 4- 
1. The principal difference with the rRBC is that the outgoing 
field at the Huygens surface location is computed using the 
first order Higdon operator. This is a better estimate than the 
teleportation operator. The outer domain is closed with a PEC 
condition. 

There are some theoretical derivations in the paper [53]. 
Especially, it is shown that the reflection produced by the ABC 
relying on the Higdon operator for estimating the outgoing 
field on the Huygens surface is identical to the reflection 
produced by the same Higdon operator used as a traditional 
ABC. The method then appears as just another implementation 
of the Higdon operators. The claimed advantage of the method 
is the easiness of stacking several surfaces to realize the 
equivalent of a high order Higdon operator. Several 
experiments show that the multiple absorbing surface permits 
an excellent absorption. However, the experiments are with a 
point source in a 2D domain where there is no evanescent 
waves, and with a scattering sphere with which results are 
provided only about the resonance frequency where the 
surrounding field is weakly evanescent. As with the rRBC, the 
experiments do not correspond to most situations of 
computational electromagnetics, because of lack of evanescent 
waves. Even by stacking a large number of absorbing surfaces, 
as long as only Higdon operators are used to estimate the 
outgoing field, the performance of the MAS method would be 
poor if the MAS surfaces were placed in the evanescent region. 


IV-3. The basis and theory of the Huygens absorbing boundary 
condition (HABC) 


The above two ABCs were developed independently and 
formulated in slightly different manners, but both rely on the 
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same basic principle that consists of cancelling the outgoing 
field leaving the domain by means of equivalent currents that 
radiate a field equal in magnitude and opposite in sign to the 
field to be cancelled (Fig. 4-1). This is a simple and attractive 
principle. These ABCs yield correct results in the test cases 
reported in the papers, but they do not significantly outperform 
such existing ABCs as the Mur or Higdon ABCs, because they 
do not absorb the evanescent fields. This prevents them to 
challenge the PML ABC in realistic calculations where there 
are always evanescent fields. 

The MAS and rRBC have been generalized in paper [56], in 
the form of an ABC called the Huygens ABC (HABC). 
Theoretical properties of the Huygens ABC are derived in [56], 
both in continuous spaces and in the discrete FDTD space. The 
theory permits the observations and results in the rRBC and 
MAS papers to be explained and interpreted, and opens the 
way to the development of new families of ABCs that may 
challenge the PML. In this section, the principle of the HABC 
and the essential theoretical facts derived in [56] are 
summarized. 

As the above MAS and rRBC, the HABC consists in 
generating an opposite wave to cancel the outgoing waves, by 
means of a Huygens surface (Fig. 4-1). In electromagnetics, a 
Huygens surface [16] permits sources to be replaced with 
currents sheets 


(4-1a) 


J5= nXH, 


RB, Sonxe. (4-1b) 
where n is the unit vector normal to the surface, oriented in the 
direction opposite to the sources, and FE; and H; are the fields 


that would exist upon the surface if the sources were present. 
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To realize a perfect HABC, E; and H; should be the opposites 
of the outgoing fields at the Huygens surface location, that is at 
the FDTD nodes where the HABC is implemented. In practice, 
things are different, because E; and H; cannot be known 
exactly. This can be viewed easily with the FDTD method. 
More fundamentally, in the continuous space this is because 
the field is discontinuous through the Huygens surface. For this 
reason, an estimate of the outgoing field is used in place of its 


exact, but unknown, value. Denoted as E (x. ,t) , the estimate 


at the Huygens surface location x, is defined as a linear 
combination of the actual field values E,(x,-6x;, t-6,) present 
at several times and several locations close to the Huygens 
surface, on its inner side. For the +x side of the Huygens 
surface, assumed as a parallelepiped box, this can be written as 


ay N 
E@a)=>. a; E_ (x, — Ox, ,t-o,) (4-2) 
k=l 
where 6x, > 0, dt, > 0, and 
N 
Sa; =1 (4-3) 


k=1 


With the FDTD method, where the field is known at mesh 
nodes separated with Ax and at times separated with At, this 
can be rewritten in terms of operators as 


E(x,,t)=PE,(x,,t) (4-4) 


with 


N 
P=X a KZ (4-5) 
k=1 

where K and Z” are the Higdon shift operators in space and 
time (see section 2-2), m, and n,; are integers (m, > 0, ng È 0, 
with mm; = ng = 0 excluded). Similar estimates as (4-2)-(4-5) are 
used for the H field, especially with the FDTD method where 
E and H are separated with Ax/2 in the direction normal to the 
Huygens surface [16]. 

The field behind the Huygens surface is the addition of the 
outgoing field with the field radiated by the Huygens surface. 


It is derived in [56] for an outgoing plane wave at incidence 0. 


It reads: 
y, Ax 
Yai mar cos 0 m, ) 
c 


k=1 


im OE incident 


behind ~~ d 


E (4-6) 


This field is reflected by the PEC ending the domain (Fig. 4-1) 
and passes backwards through the Huygens surface. Derivation 
in [56] shows that the field reflected in the inner domain reads: 


1 (4-7) 





E eftected = f E poehind dt' 


N 
Za, (n,At + cos @ m,Ax/c) 
k=l 


From (4-6) and (4-7), the Huygens surface acts as a derivator 
on time of the outgoing field, and as an integrator of the 
ingoing field propagating backwards. The waveform of the 
field reflected back into the inner domain is then the same as 
that of the outgoing wave, but its magnitude is multiplied with 
the product of the coefficients in (4-6) and (4-7). In other 


words, the HABC formed with the Huygens surface and the 
PEC ending the domain is an ABC of reflection coefficient: 


a,(cos@ m,Ax—cn,At) 


Maz 


(4-8) 


> 
un 


a,(cos@ m, Ax + cn,At) 


Maz 


> 
un 


It can be easily shown [56] that this coefficient is nothing but 
the reflection coefficient of the operator (4-5) used as a 
traditional operator ABC on the outer boundary. This was also 
derived in [53] in the special case of the Higdon operator. The 
HABC is then just another implementation of the analytical 
operator ABCs in the form (4-5). 

The teleportation operator of the rRBC is the special case of 
(4-5) Pe=K Z ' For the MAS, the estimate is computed with 
the first order Higdon operator, which reads, from (2-29) 


Pa =KZ'+wK-wZ` (4-9) 

with w from (2-20). The reflection produced by the 

teleportation operator (called the elementary operator in [56]) 

is the special case (N=1, a;=1, m;=n;=1) of (4-8): 
(cos@ Ax —c At) 


7S ey (4-10) 

(cos@ Ax +c At) 
whose zero reflection is at oblique incidence cos® = c At/Ax. In 
the rRBC [54], the outer PEC is replaced with an impedance 
condition equivalent to a first order condition, with reflection 
(2-9). The overall reflection is then the product of the two 
reflection coefficients. This permits (see [56]) the 
interpretation of some results reported in [54, 55]. For the 
MAS surface, using (4-9) into (4-8) yields, as expected, the 
reflection from a Higdon operator used as a traditional ABC, 
that is (2-9). 

Obviously, the second order Higdon operator, which is also 
a special case of (4-5), can be used with the HABC. It would 
produce a reflection equal to the square of (2-9). The 
corresponding P# operator is given in [56]. 

The results in Fig. 4-2 from [56] demonstrate the derivation 
(4-6) and re-integration (4-7) of the outgoing pulse by the 
HABC surface. The experiment was in 1D (equivalent to 3D at 
normal incidence) with the elementary operator used in the 
tRBC whose reflection is (4-10) for 6 = 0. It is clear that the 
pulse behind the HABC is the derivative of the Gaussian 
incident pulse, which is then re-integrated as it passes back 
through the HABC, with a final magnitude in accordance with 
(4-10). 

A remark should be done here about special cases of 
operator (4-5) where the first derivative is replaced with the 
second derivative in (4-6), and the integral in (4-7) replaced 
with a double integral on time. The reflection then differs from 
(4-8), but it remains equal to the reflection from the same 
operator used as a traditional ABC. More details are available 
in [56]. 

The reflection coefficient (4-8) is that of pure traveling 
waves. The reflection of evanescent waves with wavenumbers 
(2-10)-(2-11) has been derived in [56] as well. It reads 
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Fig. 4-2. Reflection from the Huygens ABC. A comparison of theory with a FDTD experiment. The wave is observed at A and B on the two sides of the HABC. 
The distances from the source to A, HABC, B, and PEC are 5, 10, 15, and 60 FDTD cells of 5 cm, respectively. 
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(4-11) 
which is again, as it can be shown, the reflection from the 
same operator used as a traditional ABC. Obviously, (4-8) is 
nothing but the special case %=0 of (4-11). In the case of 
strongly evanescent waves, where coshy >> 1, (4-11) yields 
r = 1, so that the reflection is total from the HABC with any 
operator (4-5), especially with Higdon operators. This means 
that the HABC relying on (4-5) suffers from the same 
drawback as the Higdon operator ABC, it cannot absorb the 
evanescent waves that are present in most problems of 
electromagnetics. 

Fig. 4-3 from [56] illustrates the equivalence of the HABC 
with its operator implemented as a traditional analytical ABC. 
It reports experiments with a parallel plate waveguide, using 
Higdon operators of orders 1 and 2, implemented either as an 
operator ABC or as a Huygens ABC. It can be seen that the 
results of the two implementations are superimposed, both in 
the traveling region above the cutoff frequency 3.75 GHz and 
in the evanescent region where the reflection is total (the low 
frequency plateau is 6 dB below total reflection because of the 
natural decrease of the evanescent waves over the separation 
between the source and the ABCs). 


IV-4. Implementation of the HABC in Cartesian coordinates 
and extensions of the HABC 


The Fig. 4-1 summarizes the theoretical principle of the 
HABC. However, it must be slightly modified in practice 
because the Huygens surface only radiates an estimate of the 
outgoing field. Especially in Cartesian domains where the 
HABC surface is a rectangle in 2D, or a parallelepiped in 3D. 
At the corners or edges the direction of the normal to the 
surface is discontinuous (discontinuity of orientation by 7/2). 
From this, the outgoing field estimated by means of an 
operator is discontinuous. This results in the generation of a 
spurious source of field in the vicinity of corners and edges. A 
solution has been found [57] by adding what is called 
extensions to the physical Huygens surface. They are 
represented in Fig. 4-4 in the 2D case. The four sides of the 
HABC are extended up to the outer PEC ending the domain. 
Similarly, in 3D the six faces of the HABC are extended up to 
the PEC. The Huygens currents (4-1) are applied on the 
extensions in the same way as on the normal sides of the 
HABC. With this disposal, it can be seen [57] that there is no 
discontinuity of the radiated field in the different regions 
defined by the HABC and the extensions. 

More details can be found in [57] where numerical 
experiments illustrate the need of the extensions. The practical 
implementation of the extensions is a little change in the 
encoding since is just consists of extending some loops. The 
need of extensions was not found when the work [56] was 
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done, because the experiments were only in 2D without 
corners (HABC at the end of a waveguide). In the rRBC and 
MAS papers the need of extensions is not mentioned as well. 
However, although the rationale for the authors to do this is 
not reported, it seems that in [54, 55] the Huygens planes were 
extended up to the outer boundary ending the domain in the 
numerical experiments. 
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Fig. 4-3. Comparison of the reflection coefficients computed using either an 
operator ABC or a Huygens ABC at the end of a parallel-plate waveguide. In 
both cases the operators are either the 1-order Higdon operator or the 2-order 
Higdon operator. The PEC result is with the waveguide ended with a PEC 
alone, without ABC. 
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Fig. 4-4. The extensions of the HABC. The six faces of the Huygens surface 
box (in 3D) are extended up to the outer PEC of the computational domain. 


IV-5. The potential interest of the HABC 


As shown in the above the HABC is nothing but another 
implementation of analytical ABCs, as the Mur or Higdon 
ABCs or any operator in the form (4-5). From this, one may 
think that it is of little interest. Indeed, this is not true, the 
HABC has some characteristics that render it of potential high 
interest for numerical electromagnetics. Mainly the following 
two: 

1/ As proposed in the pioneering work of the MAS [53], it 
is possible to stack several HABCs. The reflection is then the 
product of the individual reflections. Stacking for instance 
three first order Higdon HABCs is equivalent to a single third 
order Higdon operator. There is no gain, in theory, but the 
implementation is simpler, and the stability may be better, as it 
is claimed in [53]. However, the true interest in stacking 
several HABCs is not in staking numerous HABCs only 
designed for the absorption of pure traveling waves, as done in 


[53]. The true interest is indeed in the possibility of easily 
stacking operators designed for the traveling waves with 
operators designed for the evanescent waves. This is 
demonstrated in the following. 

2/ The HABC can be placed in front of another boundary 
condition, by just replacing the PEC in Fig. 4-1, with the other 
ABC. This permits easy combinations with such ABCs as the 
PML or a simple and rather exotic, but effective, real stretch 
of coordinates (see in the following). Combinations of several 
ABCs are of high potential interest because they render 
possible combinations of one ABC designed for the traveling 
waves with another designed for the evanescent waves, as 
demonstrated in the following. 


IV-6. A HABC for both traveling and evanescent waves 


In the framework of analytical ABCs, it seems only the 
Betz-Mittra operator (2-49) has been reported for the 
absorption of evanescent fields. In [56], another simple 
operator has been used to absorb the evanescent field at the 
end of a waveguide. It consists in writing that the field at the 
boundary, or at the Huygens surface node with the HABC 
implementation, is the field one cell in the interior multiplied 
with a coefficient B < 1. This is an operator designed for 
waves evanescent in the direction perpendicular to the 
boundary, with the direction of the propagation of the phase 
parallel to the boundary. In terms of operator, this reads: 


Pra =£ K (4-12) 


With the phase propagation parallel to the boundary (0 =+7/2), 
the reflection coefficient of (4-12) is derived in [56]. It reads 





® sinh y &x = 
n 1- pe" __ 1- (l+ øsinh yAx/c) (4-13) 
-2simgză 1- p(l- øsinh y Ax/c) 
1-Be< 


which is not a special case of (4-11) because (4-3) used to derive 
(4-11) does not hold with (4-12). Reflection (4-13) vanishes if 


-2 sinh Ax 

prec” 
With a parallel plate waveguide of cutoff frequency @,, the 
product @ sinhy is given by (3-21) when œ << @,. Considering 
the TM, mode and using this œ sinh% into (4-14) yields the 


value of B which cancels the reflection of the evanescent 
waves at low frequency: 


(4-14) 


pre 

Results provided in [56] show that the operator (4-12) with 

B from (4-15), implemented either as a traditional ABC or as a 
HABC, does absorb very well the evanescent waves of the 
TM, mode. The operator (4-12) can be squared to obtain a 
second-order operator (see [56]). By stacking several HABCs, 
some relying on the Higdon operator and some on (4-12), both 
the traveling waves and the evanescent waves are absorbed. 
This is illustrated in Fig. 4-5 from [56] where either two or 
four HABCs are used. The best result in the figure was 
computed with two second order HABCs for traveling waves, 
and two second order HABCs for evanescent waves, which is 
equivalent to an eight-order analytical ABC, with four orders 


Ax 


(4-15) 
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devoted to the traveling waves, and four to the evanescent 
waves. The reflection is then roughly similar to that obtained 
with an eight-cell optimized CFS-PML. 
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Fig. 4-5. Reflection from the end of a waveguide using various combinations 
of Higdon operators (first order Phigl and second order Phig2) with 
evanescent operators (4-12) (first order Peval and second order Peva2). 


IV-7. Combination of the HABC with the PML 


It is easy to place a HABC in front of a PML ABC. An 
experiment is reported in [56] with the same waveguide as 
above. Here the HABC is used to absorb the TMy mode of the 
parallel plates which is a traveling mode at any frequency (no 
cutoff). The CFS-PML optimized for the TM; mode cannot 
absorb the traveling waves below the cutoff (it is there just a 
real stretch of coordinate). In principle, this situation could be 
addressed by adding some cells of PML with a normal 
stretching factor to absorb the TMy mode. However, a more 
elegant solution consists in adding a HABC with the Higdon 
operator in front of the PML. This is simpler and probably less 
demanding in computational resource. The reader is referred 
to [56] to see more details on this experiment. 


IV-8. The stretched-mesh HABC (SM-HABC) 


The SM-HABC [58] is another illustration of the possibility 
of combining the HABC with another ABC. It consists in 
using a HABC for the absorption of traveling waves and a 
stretched mesh (SM) for the absorption of evanescent waves. 
The SM ABC is not a new ABC, it is just a strongly stretched 
FDTD mesh used to fill a large enough domain outside the 
HABC, in such a way that the evanescent waves can decrease 
to a negligible level. Stretched meshes were used in the past to 
reduce the overall number of FDTD cells while preserving a 
large physical domain. However, the stretching was severely 
limited when using ABCs on the outer boundary, because of 
the reflection of high frequencies from the large cells. 
Conversely, by placing the HABC in front of the stretched 
mesh, the high frequency traveling waves are absorbed by the 
HABC before entering the stretched mesh (more precisely, the 
high frequencies are reflected by the stretched mesh, which 
plays the role of the outer PEC in Fig. 4-1). Only evanescent 
waves are present behind the HABC, in the stretched domain. 
Since in most cases the evanescent waves are low frequency 
waves, a strong stretch can be applied without penalty, so that 


the stretched mesh can be reduced to a few FDTD cells while 
preserving a large enough physical domain. The principle of 
the SM-HABC is represented in the Fig.4-6. A scattering 
object is surrounded with the HABC and its extensions. A 
coarse mesh fills the space outside the closed part of the 
HABC. 

The SM-HABC is well suited to such problems as the 
interaction of an incident wave with a scattering structure, 
where the high frequencies are traveling waves and the low 
frequencies are evanescent waves, with a transition about the 
fundamental resonance of the structure. The high frequencies 
are then absorbed by the HABC which relies on the first order 
Higdon operator in [58], and the low frequencies decrease in a 
natural manner in the large domain outside the HABC. 
Experiments with this kind of problems are reported in [58]. 
However, the concept could be more widely applied in 
numerical electromagnetics since in most problems the 
evanescent waves are low frequency waves, allowing then 
their attenuation by means of a stretched mesh. 


i 
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HABC 

n 
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mesh un mesh . 
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Fig. 4-6. The principle of the SM-HABC. The HABC (closed surface) is 
surrounded with a stretched mesh (geometrical expansion of the FDTD cell). 
The HABC is placed dy4gc from the object and the outer PEC dpzc from the 
object. 


The experiments reported in [58] shows that the Higdon- 
based HABC can very well absorb the traveling waves, even 
when it is placed in the close vicinity of a scattering structure. 
This is illustrated by the Fig. 4-7 where the HABC is placed 
several distances from a 100-10-0-cell PEC plate stricken by 
an incident pulse (a Gaussian pulse or a Unit-step pulse). A 
large domain filled with a uniform mesh surrounds the HABC 
in order that the evanescent waves decrease to a negligible 
value (dpgc = 200 cells is the distance between the 100-10-0- 
cell plate and the outer PEC). Even with the HABC placed 
only three FDTD cells away the object (daagc = 3 cells in the 
figure), there is no visible reflection of the traveling waves 
(taking account of the delay of propagation over 400 cells of 1 
cm, if a reflection were present it would appear about time 20 
ns in the figure). These results show that the first order Higdon 
ABC can achieve a very good attenuation of traveling waves, 
better than was previously believed. The HABC 
implementation permits the absorptions of traveling and 
evanescent waves to be decoupled. This clearly shows the 
high performance of the operator, conversely to the case 
where it is implemented as a traditional ABC, on the outer 
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boundary, where the reflections of traveling and evanescent 
waves cannot be clearly distinguished. 

In [58], the cell size in the stretched mesh grows 
geometrically from the HABC to the outer PEC, i.e. the cell is 
multiplied with a factor g from one cell to the next. The stretch 
can be large, resulting in a dramatic reduction of the number 
of cells needed to fill the domain outside the HABC. This is 
illustrated by the Fig. 4-8 where the same problem as in Fig. 
4-7 is addressed, i.e. the field on the surface of a 100-10-0-cell 
plate. The HABC is three cells from the plate, and ng stretched 
cells outside the HABC replace the uniform mesh of Fig. 4-7 
(the physical domain is identical, equivalent to 200 uniform 
cells). As can be seen, the results remain acceptable, even with 
a small number of stretched cells in place of the large number 
of uniform cells. In other words, 3 strongly stretched cells can 
replace 196 uniform cells (the stretch begins one cell away 
from the HABC which is 3 cells from the object). The 
computational requirements, memory and CPU times, are then 
dramatically reduced in comparison with the uniform mesh. 
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Fig. 4-7. Electric field on the surface of a 100-10-0-cell plate, computed with 
the HABC placed 3, 10, 30, 100 cells from the plate, a uniform mesh in the 
whole domain, and the outer PEC 200 cells from the plate. 
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Fig. 4-8. Electric field on a 100-10-0-cell plate computed with the HABC 3 
cells from the plate (dyagc = 3 cells) and four stretched meshes (ng = 3, 5, 7, 


10 cells) behind the HABC (mesh 2 defined in [58]). In all the cases the 
plate-PEC separation equals 2 m which is equivalent to 200 uniform cells. 


The achievement in Fig. 4-8 is possible because the high 
frequency traveling waves are absorbed by the Higdon-based 


HABC in front of the stretched mesh. With the traditional 
implementation as an analytical ABC on the outer boundary, 
using a widely stretched mesh is not possible, since then the 
high frequency waves cannot propagate through the large 
FDTD cells, and then cannot reach the Higdon operator. This 
is clear in the Fig. 4-9, where the same problem as in Fig. 4-8 
is addressed, except that the HABC is replaced with the 
Higdon ABC implemented on the outer boundary, in place of 
the outer PEC. Then the high frequencies are reflected 
backwards by the stretched mesh, resulting in the strong 
spurious reflection visible in Fig. 4-9. Thus, comparing Fig. 4- 
8 with Fig. 4-9 clearly demonstrates the high interest of the 
HABC implementation of operator ABCs in comparison with 
the traditional implementation. The possibility of placing the 
HABC in front of another ABC is a unique feature that opens 
the way to new achievements in the matter of ABCs. 

An example of application of the SM-HABC to the airplane 
problem previously addressed with the Higdon operator in Fig. 
2-7 and with the CFS-PML in Fig. 3-10 is reproduced in Fig. 
4-10. It can be seen that the computed results are correct with 
a small number of stretched cells (ng = 3-5). The overall 
computational domain size and CPU time are very similar to 
the ones of the CFS-PML calculations in Fig. 3-10. The SM- 
HABC cannot really outperform the CFS-PML, but it can 
challenge it, with the advantage of a simpler implementation. 


— Reference 
— — ng=3 cells 
Mesh 2 H dHigdon <=> 200 cells H 777 "9= 5cels | 
ARAR ng = 7 cells 
ng = 10 cells | 








2000 \ 



















1500 


1000 










500 Wat We eats tie 
i Mahi Wital nf 
Eito jip! hy t My ps 
hatina bi gh Adva Alil 
A i bis «| 








E Field normalized to Incident Field 


-500 











0 5 10 15 20 25 30 35 40 45 50 55 60 
Time (ns) 


Fig. 4-9. Electric field on a 100-10-0-cell plate computed with four stretched 
meshes (the same as in Fig. 4-8) and a Higdon operator ABC placed on the 
outer boundary. In all the cases the plate-ABC separation equals 2 m which is 
equivalent to 200 uniform cells. 


IV-9. Summary and discussion on the HABC 


Although it is in principle just another implementation of 
analytical ABCs, the HABC is much more for two reasons. 
Firstly, it allows easy combinations of operators devoted to 
traveling waves with operators devoted to evanescent waves. 
Secondly, it can be placed in front of another ABC and then 
combined with a PML, or with more innovative ABCs such as 
the simple, but effective, real stretch of coordinates. 

From this, novel ABCs that can really challenge the PML 
ABC can be elaborated. As the SM-HABC introduced in [58] 
for wave-structure interaction problems, but that could be 
applied as well to other problems where evanescent waves are 
present at low frequency (as waveguide problems for 
instance). Moreover, it seems that the SM-HABC has rather 
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more potential than the CFS-PML. Because the optimization 
of the CFS-PML is more problem dependant. It assumes that 
the frequency spectrum is like in Fig. 3-9, with the optimum 
parameter ©, depending on the relationship (3-21) for 
waveguide problems, or (3-24) for interaction problems. In 
both cases, the relationship depends on a parameter, the mode 
in (3-21) or the parameter p in (3-24). Conversely the SM- 
HABC only relies on the fact that evanescent fields are low 
frequency fields, so that it can be suited to all the modes in a 
waveguide problem and can even absorb low frequency 
traveling waves (as the TM, mode). Although the CFS-PML 
could in principle address any situation by varying its O% 
parameter within it, the SM-HABC appears as superior. More 
numerical experiments should be performed to conclude and 
compare more extensively the two ABCs in situations close to 
realistic applications, but the less restrictive assumptions of 
the SM-HABC are, a priori, a favorable feature. 
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Fig. 4-10. Electric field on the surface of a 504-632-160-cell airplane with 
SM-HABCs of 3, 5, 7 stretched cells, for a unit-step and a double-exponential 
pulse. 
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On the other side, the HABC is only described in [56-58] 
for the simulation of free space in a vacuum. Obviously it 
could also be used with little change in isotropic and non 
dispersive media. However, its extension to more complex 
media is far from trivial. Finding analytical ABCs for general 
media has been a challenging and mainly unsolved question in 
the past. Thus, at present there is no available operator to 
estimate the outgoing field in complex media. This is a serious 


limitation of the HABC in contrast with the PML which has 
been relatively easily extended to any kind of media, from 
lossy media to metamaterials. 


V. CONCLUSIONS AND PROSPECTS 


Absorbing boundary conditions are needed to solve most 
problems of numerical electromagnetics with such finite 
methods as the FDTD or the FEM. They have been the subject 
of researches almost continuously over the past 40-45 years. 
Several analytical ABCs were developed in the years 1970- 
1990. However, despite numerous works, little progress was 
observed during this period. The PML ABC introduced in 
1994 rapidly appeared as a breakthrough advance because it 
can outperform previous ABCs in the two principal aspects of 
the simulation of free space: 1/ the reduction of the reflection 
from the ABC, which means a better accuracy of the 
calculations, and 2/ the reduction of the overall computational 
cost of the calculations. In addition, the PML has been 
extended easily to any media, including dispersive, 
anisotropic, and metamaterial media. However, even if the 
PML is effective and universal, the more recently introduced 
Huygens ABC enriches the range of available ABCs. It 
permits re-use of the existing analytical ABCs in a different 
and successful manner, and subsequently permits the design of 
ABCs that can challenge the PML, and may outperform it in 
some problems, especially regarding the question of the 
simplicity. 

What is clear nowadays is that the problem of the 
simulation of free space has not been well understood in the 
years 70 and 80. The researches were focused on the 
improvement of the absorption of traveling waves, as if only 
traveling wave existed. In reality, evanescent waves are 
present in most problems. This explains why there has not 
been significant progress during this period. The reason of the 
effectiveness of the PML and HABC mainly lies in their 
ability to absorb the evanescent fields. With the PML by 
means of complementary reflections (Ramahi effect), or by 
using an optimized CFS stretching factor. With the HABC by 
means of operators designed to this end, or by means of 
combinations with such non conventional ABCs as a simple 
real stretch of coordinates. 

What may explain why the workers have not paid much 
attention to the evanescent waves in the past decades is the 
absence of the general evanescent solutions of the Maxwell 
equations in the textbooks on electromagnetism, for instance 
in Stratton or Van Bladel textbooks. This lack of the general 
waveform of the electromagnetic field has been an obstacle to 
interpret the reflection from PMLs, in the context of which the 
evanescent solutions of the Maxwell equations had to be 
derived from scratch [40] in the 2D case (appendix A). This 
probably has been an obstacle as well for the development of 
effective analytical ABCs in the 70s and 80s. 

Finally, although nowadays the simulation of free space is 
considered as satisfactory in most situations, mainly using the 
PML ABC, progress is endless and there is still room for 
works in the field of ABCs. With such objectives as the 
simplicity of the implementation, the simplicity of the settings 
of the ABC, or the reduction of the sensitivity of the results to 
the parameters of the problem. In a vacuum, which is the most 
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current medium in applications, the Huygens ABC opens the 
way to new researches, mainly because of its unique feature 
that permits easy combinations of different ABCs. However, 
further researches can only succeed if one keeps in mind the 
central question in the matter of ABC research, which is, once 
again, the fact that the evanescent waves have to be absorbed. 
This in turn implies that a clear and explicit knowledge of the 
general solutions of the Maxwell equations is available. 


APPENDIX A 


EVANESCENT SOLUTIONS OF THE MAXWELL EQUATIONS 


Let us consider the equation of dispersion of the Maxwell 


equations in the 2D case: 
2 


kè+ k=? (A-1) 





2 


Beside the usual solution (k, = @cos@/c, k, = œ sin0/c), the 
following more general solution also satisfies the equation 


k, = (cosh y cos0+ jsinh 7 sin 6) (A-2) 
c 


k, = (cosh ysin 0- jsinh ycos6) (A-3) 
where 0 and x are free parameters. The waveform exp (jat- 


Jk,x- jkyy) of the E (or H) field is then in the form 


_ cosh y 


E=E, et c 





(cos 4+ yin 2) | -2 sinh y (y cos 0- xsin 0) 
es 


(A-4) 


By the rotation of coordinates defined in the figure, this can be 
rewritten as 


cosh 
-eZ y | -2 sinh zY 
: 


poe ‘ e 


(A-5) 





Fig. A-1. Rotation of coordinates 0. 


Equation (A-5) is the waveform of a plane wave whose 
phase propagates with speed cohy/c in direction X forming the 
angle © with direction x, and whose magnitude decreases in 
direction Y perpendicular to X. Quantity yx € [-00, +00] is the 
evanescence parameter. The case x = 0 is the special case of 
the pure traveling waves. 

Some results are derived in [40] about the reflection of the 
evanescent waves from an interface, which also hold for the 
reflection by an ABC. For an incident wave whose phase 
propagates in direction @ the phase of the reflected wave 
propagates in direction 7-9, as that of a pure traveling wave. 
And the evanescence parameter x of the reflected wave is 
equal in magnitude, but opposite in sign, to that of the incident 
wave. 


The 2D solutions of the Maxwell equations are sufficient to 
compute the reflection coefficients from an analytical ABC or 
from a PML, because ABC problems are mainly 2D problems 
(at least as long as the medium is isotropic). 


APPENDIX B 


EVANESCENT FIELDS IN THE VICINITY OF SCATTERING 
STRUCTURES 


Apparently, in the literature or textbooks, there is nothing 
about the general form of the evanescent fields surrounding 
structures stricken by an incident wave. Only canonical cases 
are addressed. In order to interpret the observed reflection 
from PMLs, it was assumed in the work [40] that the low 
frequency field around PEC structures is composed of 
evanescent waves whose evanescence parameter % is 
connected with the frequency by the relationship: 


@ sinh y= p— (B-1) 

w 
where p is a parameter of the order of unity and w is the 
(largest) size of the structure. The relationship was deduced 
from the fact that the evanescent fields strongly decrease over 
a distance of the order of the largest size of scattering 
structures. Taking account of the coefficient of decrease of the 
evanescent waves (appendix A, formula A-5), this reasoning 
yields: 


oO. 
——sinh y w 
c 


e =e? (B-2) 


or equivalently (B-1). The formula (B-1) is fundamental to 
optimize the CFS-PML (its Œy parameter). It could also be 
used to design operators for the absorption of evanescent 
waves in scattering calculations. 

Formula (B-1) can be viewed as obtained empirically. 
However, more recently the question has been investigated 
theoretically and a rigorous formula derived in the special case 
of symmetric 2D structures [59]. In that case, the spectrum of 
the plane waves evanescent in the direction perpendicular to 
the largest size w presents peak values for the following 
% parameters 


o cosh y, = 2m+1) 2E (B-3) 
w 


with magnitudes of the peaks decreasing as m grows. Since the 
resonance frequency is ©% = T c/w, far below the resonance 
cosh% >> 1 holds. From which cosh% ~ sinh% and then (B-3) 
agrees well with (B-1). It can be noticed that (B-3) is the same 
as the formula that holds for a parallel plate, with just the 
object size w in place of the waveguide width. 
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